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ABSTRACT 


In  many  applications  involving  wave  propagation,  problem  domains  are  often 
very  large  or  unbounded.  A  common  numerical  method  used  to  solve  such  problems 
is  to  truncate  the  domain  via  artificial  boundaries  to  form  a  finite  computational 
domain.  To  accomplish  this,  Non- Reflecting  Boundary  Conditions  (NRBC’s)  which 
minimize  spurious  wave  reflections  are  imposed,  The  quality  of  the  solution  strongly 
depends  on  the  properties  of  both  the  NRBC  and  the  wave  behavior. 

This  dissertation  explores  the  use  of  Higdon  NRBC’s  to  solve  shallow  water 
equations  (SWE’s)  in  a  dispersive  environment.  A  linearized  SWE  model  is  developed 
that  includes  stratification  and  advection  effects.  Initially  a  single  NRBC  is  used  to 
truncate  a  semi-infinite  channel.  Later  four  NRBC’s  are  used  to  restrict  an  infinite 
plane.  In  both  cases  finite  rectangular  domains  are  formed.  A  scheme  developed  by 
N eta  and  Givoli  is  used  to  rapidly  discretize  high- order  Higdon  NRBC’s.  Finite  differ¬ 
ence  methods  and  are  used  in  all  numeric al  schemes,  which  are  solved  explicitly  when 
possible.  Results  will  show  that  Higdon  NRBC’s  can  be  used  effectively  to  restrict 
large  rectangular  domains  when  solving  SWE’s  that  include  the  before  mentioned 
effects. 
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I.  INTRODUCTION 

Phenomena  involving  the  propagation  of  waves  in  unbounded  (or  very  large) 
domains  are  applicable  to  many  fields  including  acoustics,  electromagnetics,  meteo¬ 
rology,  and  geophysics.  However,  it  is  infeasible  to  compute  numerical  solutions  for 
regions  of  this  scope.  Therefore,  it  is  necessary  to  define  artificial  boundaries  that 
reduce  the  size  of  the  domain.  To  accurately  model  the  wave  action  in  the  truncated 
region,  artificial  boundary  conditions  must  be  imposed  that  allow  waves  propagating 
inside  the  region  to  pass  freely  without  spurious  reflections,  which  would  pollute  the 
computational  domain.  Such  a  boundary  condition  is  known  as  the  Non- Reflecting 
Boundary  Condition  (NRBC)  and  is  the  main  subject  of  this  dissertation. 

In  general,  it  is  not  possible  to  construct  a  boundary  condition  that  will  accom¬ 
plish  the  criteria  perfectly,  but  during  the  last  25  years  research  has  been  conducted 
to  develop  NRBC’s  that  after  discretization  lead  to  stable,  accurate,  efficient  and 
easily-implemented  schemes  [Ref.  1,  2,  3,  4],  Investigations  in  the  late  70’s  to  early 
30’s  produced  a  number  of  low-order  local  NRBC’s,  e.g,  the  Engquist-Majda  [Ref. 
5]  and  Bayliss-Turkel[Ref,  6]  boundary  conditions.  The  exact  non-local  Dirichlet-to- 
Neumann  (DtN)  NRBC  [Ref.  7,  8]  and  the  Perfectly  Matched  Layer  [Ref.  9]  boundary 
conditions  were  developed  in  the  late  30’s  and  early  90’ s.  Subsequently,  higher  order 
NRBC’s  were  introduced,  but  were  difficult  to  employ  beyond  the  2nd  or  3rd  order. 

Only  since  the  mid  90’s  have  practical  higher  order  schemes  been  developed. 
Collino  [Ref.  10]  proposed  such  a  scheme  for  two-dimensional  time- dependent  wave 
in  a  rectangular  domain.  Grote  and  Keller  [Ref.  11]  extended  the  domain  to  three 
dimensions  in  a  scheme  based  on  spherical  harmonic  transformations.  They  extended 
their  work  to  include  elastic  waves  [Ref.  12],  These  findings  were  independently  pub¬ 
lished  by  Sofronov  [Ref.  13]  in  Russian  literature.  Hagstrom  and  Hariharan  [Ref.  14] 
constructed  high-order  NRBC’s  for  two-  and  three- dimensional  domains  based  on  the 
analytic  series  representation  for  the  outgoing  solutions  of  these  equations.  Guddati 
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ajid  Tassoulas  [FLef .  15]  devised  a  high- order  NRBC  scheme  foi  time- dependent  waves 
in  a  2- dimensional  wave  guide  using  rational  approximation  and  recursive  continued 
fractions.  Givoli  [FLef .  16]  derived  high-order  NRBC’s  for  ageneral  class  of  wave  prob¬ 
lems  leading  to  a  symmetric  finite  element  formulation.  These  early  investigations 
utilized  either  time-harmonic  waves  or  non- dispersive  time- depen  dent  waves. 

Wave  dispersion,  however,  is  an  ever  present  phenomenon.  In  the  late  80 ’s  and 
early  90’s,  Higdon  developed  NRBC’s  for  non- dispersive  waves  [Ref.  17,  18,  19,  20], 
but  later  showed  that  his  schemes  could  be  applied  to  the  dispersive  (Klein- Gordon) 
wave  equation  [Ref.  21].  Higdon’s  work  involves  low  order  formulation  of  his  scheme. 
Givoli  and  Neta  [Ref.  22,  23,  24]  present  an  algorithm  for  implementing  the  Higdon 
NRBC  to  any  order  using  high- order  FD  discretization.  They  further  developed 
methods  to  rewrite  the  Higdon  NRBC  without  using  high  order  derivatives  and  to 
generate  Higdon  parameters  that  maximize  the  non-reflection  property  of  the  NRBC 
in  a  dispersive  wave  environment. 

In  the  present  work  I  will  develop  high  order  Higdon  NRBC  schemes  for  use 
with  linearized  shallow  water  equations  (SWE’s)  in  Cartesian  coordinates.  The  SWE 
model  is  further  enhanced  to  include  the  effects  of  stratification  and  advection,  A 
single  Higdon  NRBC  is  initially  applied  as  an  artificial  boundary  on  one  side  of 
a  semi- infinite  channel.  Later  four  Higdon  NRBC’s  are  applied  to  the  sides  of  a 
rectangular  domain  to  restrict  an  infinite  plane.  Finite- difference  schemes  are  used  to 
numerically  solve  the  problems.  Discrete  forms  of  the  Higdon  NRBC,  based  on  the 
work  of  Givoli  and  Neta,  are  then  employed  on  the  artificial  boundary,  The  results  of 
several  numerical  examples  are  reported  to  validate  the  SWE  models  as  well  as  the 
use  of  the  Higdon  NRBC  as  an  effective  means  of  restricting  a  very  large  domain. 
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II.  MODELING  GEOPHYSICAL  FLUID 

FLOW 

In  this  chapter  we  develop  a  model  that  describes  large  scale  geophysical  flow. 
Many  of  the  details  of  this  derivation  are  developed  by  Pedlosky  [Ref.  25],  We 
start  by  considering  the  dynamics  of  a  shallow  rotating  fluid  layer.  The  fundamental 
condition  that  characterizes  a  shallow  layer  is: 


where  D  and  L  characterize  the  scale  of  vertical  and  horizontal  motions.  This  charac¬ 
terization  is  applicable  to  large  scale  atmospheric  and  oceanic  flows  where  the  vertical 
scale  of  the  fluid  layer  is  of  the  order  of  miles,  while  the  horizontal  scale  is  of  the  order 
of  hundreds  or  thousands  of  miles.  Work  by  Rossbv  [Ref.  30]  and  Stommel  [Ref.  31] 
show  that  such  a  working  geophysical  model  can  be  developed  by  assuming  that  the 
fluid  is: 

•  Incompressible  (density  independent  of  pressure), 

•  Inviscid  (no  internal  frictional  forces), 

•  Homogeneous  (not  stratified  with  regards  to  density). 

Later  in  this  paper,  the  “homogenous75  assumption  is  relaxed.  The  following  physical 
laws  are  applied  in  deriving  a  geophysical  model: 

•  Conservation  of  Mass, 

•  Conservation  of  Momentum  (Newton’s  2f3d  Law), 

•  Conservation  of  Energy, 

•  Second  Law  of  Thermodynamics. 

Since  the  fluid  is  incompressible,  the  energy  equations  are  uncoupled  from  the  model. 
Hence,  we  consider  only  the  conservation  of  mass  and  momentum,  The  basic  form 
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□f  these  laws  apply  to  fixed  quantities  of  mattei.  In  the  analysis  of  fluid  flow,  we 
are  concerned  with  fixed  volumes  (e.g.  the  properties  of  a  fluid  within  a  control 
volume  rather  than  properties  of  individual  particles  in  motion).  What  follows  is 
a  derivation  of  the  conservation  of  mass  and  momentum  equations  as  they  apply  to 
control  volumes. 

A.  CONSERVATION  OF  MASS  EQUATIONS  FOR 
FLUIDS  IN  A  CONTROL  VOLUME 

The  conservation  of  mass  law  states  that  the  total  change  of  mass  in  a  control 
volume  must  be  equal  to  the  net  flow  of  mass  entering  and  leaving  at  the  volume 
surface.  If  p  is  the  density  of  the  fluid  and  pu  ■  n  is  the  mass  flux  at  a  point  on  the 
volume  surface,  then  by  invoking  the  conservation  of  mass  law  we  can  write: 

d  f  f 

j  PdV  =  -  j(pu-n)dA,  (III) 

where  V’  is  the  volume  and  A  is  the  volume  surface  area.  By  the  divergence  theorem 
we  know  that: 

j  (pu-  n  )dA  =  j  V  •  (pu)dV.  (II. 2) 

Applying  this  to  (II.l)  yields: 

/  (|£  +  V  •  (/»))  dV  =  0,  (II  .3) 

which  implies: 

^+V-(pi)=D.  (II.4) 

Since  the  fluid  is  assumed  to  be  incompressible  and  homogeneous  (e.g.  p  is  constant), 
this  equation  is  rewritten  as: 

V  ■«=£),  (II. 5) 

or,  in  three  dimensional  Cartesian  coordinates: 

(II 


dv.  dv  dw 

ai  +  ai+  al=0’ 
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where  v,  r,  and  w  are  the  velocity  components  in  the  x-,  y-,  and  ^-directions  respec¬ 
tively.  This  expression  is  called  the  continuity  equation.  It  describes  the  conservation 
of  mass  for  a  volume  element  containing  an  incompressible  fluid  and  is  the  first  con¬ 
trolling  equation  of  the  geophysical  fluid  flow  model. 

B.  THE  MOMENTUM  EQUATIONS  FOR  FLUIDS  IN  A 
CONTROL  VOLUME 

The  next  set  of  controlling  equations  is  the  momentum  equation  and  is  based 
on  Newton’s  Second  Law: 


Force  =  ( mass)(acceleration ). 


Acceleration  components  must  be  carefully  derived,  because  geophysical  phen omen om 
occur  in  a  non- inertial  rotating  frame.  In  the  following  analysis,  we  first  consider  a 
derivation  of  these  equations  for  a  fluid  in  an  inertial  two-dimensional  Cartesian  sys¬ 
tem.  A  rotational  component  is  then  introduced  to  generate  non- inertial  momentum 
terms.  A  third  dimension  is  added  and  acceleration  components  in  vector  form  are 
obtained.  Finally  an  Earth  model  is  developed  on  which  the  geophysical  fluid  flow 
equations  of  momentum  will  be  based. 


1.  Acceleration.  Components  for  a  Fluid  in  a  Cartesian 
Inertial  Frame 

We  first  consider  acceleration  in  an  inertial  frame  (e.g.  one  that  is  not  accel¬ 
erating  or  rotating).  The  acceleration  component  in  the  x- direction  is: 

dv 

°*=&’  (IU) 


where  v  is  the  velocity  in  the  x-direction.  Since  we  do  not  treat  fluids  as  individual 
particles,  but  rather  as  a  continuum  of  matter,  v  not  only  depends  on  time,  but  also 
on  the  spatial  components  x,  y,  and  Therefore: 


dv  dv  dv  dv 
dv.  —  dt  -\-  d  x  -\-  dy  -(-  > 


(II  S) 
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ajid  (II. 7)  becomes: 


du  du  dudx  du  dy  du  dz 

**  dt  dt  ^  dx  dt  ^  dy  dt  ^  dz  dt’ 


which  is  .rewritten  as: 


du  du  du 


du 


(II.9) 


(1110) 


where  v  and  w  are  components  of  velocity  in  the  y-  and  z- direction  respectively.  We 
now  define  the  following  special  symbol: 


D  d  d  d  d 

«  =  a  +  *&  + 1 


(ii.il) 


which  is  often  called  the  derivative  following  a  fluid  [Ref.  27],  Using  this  symbol  we 
write  acceleration  as: 


a*  = 


Du 

Dt' 


Similarly: 


or  in  vector  form: 


Dv  Dm 

<lv  =  dJld  dj  ■—  —  . 

v  Dt  Dt’ 


a  = 


Du 

~Dt' 


(11.12) 


(11.13) 


(11.14) 


We  now  consider  the  effects  of  the  rotation  on  the  acceleration  vector. 


2.  Acceleration  Components  in  a  Two-Dimension 
Cartesian  Rotating  Frame 

Rotation  is  an  important  factor  in  dynamics  if  the  time  iu.  to  complete  a 

single  rotation  is  on  the  order  of  or  less  than  the  time  tL  taken  by  an  object/fluid 

/  L\ 

field  to  cover  a  distance  L  at  a  speed  U  ( e.g,  —  <  1  where  tL  =  —  J .  For  example, 

consider  a  large  disc  that  completes  one  rotation  every  5,000  seconds  (i^,.).  On  this 

disc  a  particle  travels  10  kilometers  (L)  at  a  speed  of  1  meters-second-1  (Z7),  it  will 

complete  its  journey  101  seconds  (tz).  Thus  —  =  .5  and  ■we  conclude  that  rotation 

*z 

will  influence  the  momentum  equations. 


S 


4 


Figure  1,  Cartesian  Two-Dimensional  Rotating  Frame 

In  deriving  the  rotational  acceleration  terms  for  the  simplified  two-dimensional 
model,  tracking  variables  is  a  challenging  task,  Figure  1,  which  superimposes  a  two- 
dimensional  Cartesian  rotational  frame  on  an  inertial  frame  at  a  common  origin, 
introduces  these  variables.  Additionally,  Table  I  tabulates  the  variables  with  regards 
to  their  reference  frames, 

At  time  t,  the  rotating  j-axis  makes  an  angle  fit  with  the  inertial  X-axis. 


Table  I.  Variables  Used  to  Derive  Rotational  Components  of  Momentum  Equation 


Inertial  Frame 

Rotational  Frame 

ft 

angular  rate  of  rotation 

Unit  Vectors 

I,  J 

i,  j 

Coordinates 

-*,  Y 

y 

Velocity  Vector 

U 

u 

Velocity  Components 

U,  V 

u,  r 

Acceleration  Vector 

A 

a 

Acceleration  Components 

A,  B 

a,  b 
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Applying  equations  of  rotation,  it  follows  that: 


i  =  I  cos  Q.t  -1-  J  sin  fit , 
j  =  —  I  sin  fit  -|-  J  cos  fit, 

and: 

x  =  X  cos  fit  -f  Y  sin  fit, 

y  =  —X  sin  fit  -f  Y  cos  fit. 

The  velocity  vector  in  the  rotational  frame  is: 

dxT  dyr  -  r 
u=*,+  *J=”+lJ' 


(11.15) 


(11.16) 


(11.17) 


where  the  first  time  derivative  yields  the  velocity  components  of  the  rotational  frame: 


dx  dX  dY 

u  =  —  =  - cos  fit -I- - sin  fit  —  fi-Y  sin  fit -|- fiT  cos  fit, 

dt  dt  dt 


dv  dX  dY 

v  =  —  =  —  — —  sin  fit  -f  - —  cos  fit  —  fi.Y  cos  fit  —  flT  sin  fit. 
dt  dt  dt 


(11.18) 


Similarly,  the  velocity  vector  in  the  inertial  frame  can  be  expressed  as: 


dX~  dY  - 

u=  dTI+ irJ  =  t/I+v'J’ 


(11,19) 


Algebraic  and  trigonometric  manipulation  of  (11.15)  yield: 


I  =  i  cos  fit  —  j  sin  fit, 
J  =  isin  fit  -|-  j  cos  fit. 


Using  these  and  (11.19)  yields: 
U  = 


cos  fit  4-  sin  fit  i  -f 
dt  dt 


dX  .  dY  ) 

■— —  sm  fit  -f  — —  cosfit 
dt  dt  j 


(11.20) 


(11.21) 


This  expression  along  with  (11.15)  and  (11.18)  reveal  the  following  relationships  be¬ 
tween  the  inertial  and  rotational  velocities: 


U  =  v.  —  fly, 
V  =  v  -f  fix. 


(11,22) 
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We  now  consider  the  rotational  components  of  acceleration: 


d2xr  d2yr 


a  = 


dt*'1  dt2' 

Taking  the  second  derivative  with  respect  to  time  of  (11,16)  yields: 


a  = 


!d2X  ^  d2Y  . 

cos  fit  4-  —~r  sin  fit 
dt2  dt1 


4  2f 1 
”|2 


dX  dY 

-——sin  fit  4  — —  cos  fit 
dt  dt 


—  fl2  (—X  cos  fit  4  fly  sin  fit) , 


6  = 


d2X  .  ^  d2Y  \  fdX  ^  dY  . 

sin  fit  4  “tv  cos  fl*  ]  —  2fl  — —  cos  fit  4  — —  sm  fit 

f  \  dt  dt 


dt2  dt2 

^  ^  t 

—  fl2  (— X  sin  fit  4  flV  cos  fit) . 
We  express  the  inertial  acceleration  vector  as: 


d2Y-  d2Y  ~ 

A  =  ^I+^J  =  AI+BJ' 


This  and  (11,20)  result  in: 


-  (d2X  d2Y  .  \  -  /  d2X  .  _  d2Y 

A  =  [  —  cos  fit  H — — —  sin  fit  i  4 


dt2 


dt2 


dt2 


sm  Sit  H — —  cos  fit  l . 
dt2 


(11.23) 


(11.24) 


(11,25) 


(11,26) 


Based  on  the  preceeding,  we  generate  a  final  simplifying  expression  that  relates  the 
rotational  and  inertial  components  of  acceleration: 


a  =  A  4  2flV"  —  fl2x, 
b  =  B  4  2 fit/  —  fl2y, 

or  combining  this  with  (11,22)  we  have; 


(11,27) 


A  =  a  —  2  fir  —  fl2x, 
B  =  64  2fltr  —  fl2?/, 


(11,28) 


This  equation  shows  the  effects  of  rotation  on  the  inertial  acceleration  vector.  The  first 
contribution  proportional  to  fl  and  velodfy  is  called  Coriolis  acceleration.  The  second 
contribution  proportional  to  fl2  and  the  coordinates  is  called  centrifugal  acceleration. 
We  now  refine  these  equations  for  the  Earth  model. 
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3.  Acceleration  Components  in  the  Earth  Model 

From  henceforth  we  shah  use  the  terms  geophysical  and  “large-scale”  inter¬ 
changeably,  but  to  be  more  precise  a  quantitative  definition  is  offered.  Let: 

_ _ 2tt _ 

time  to  complete  one  revolution/rotation " 

For  example,  =  7.27*10-5  radians-second-1  (one  rotation  every  24  hours).  As 

before,  we  define: 


L  =  the  length  of  a  phenomenom, 

V  =  the  speed  of  a  particle  within  the  phenomenom. 

Using  L  and  U  we  introduce  a  dimensionless  parameter  called  the  Rossby  number: 


U 

onL’ 


(11.29) 


If  the  Rnssby  number  €  is  of  order  one  or  less,  the  phenomenom  is  barge-scale”  and 
the  Earth’s  rotation  is  a  significant  factor  in  the  momentum  equations  [Ref.  25],  For 
example,  a  500  kilometer  long  ocean  current  with  a  speed  of  15  meter-sec  on  d-L  has 
a  Rossby  number  €  =  .205.  Therefore  the  current  is  large-scale  (or  geophysical)  and 
is  significantly  affected  by  Earth’s  rotation. 

In  order  to  model  the  Earth,  the  momentum  equations  for  a  rotating  sphere 
are  now  considered.  This  non-in ertdal  frame  is  complicated  by  the  fact  that  at  any 
given  point  on  the  sphere  we  perceive  ourselves  on  a  planar  disc.  To  recreate  this 
“human  experience”  we  set  up  a  three  dimensional  Cartesian  system  whose  origin  is 
at  the  point  of  the  observer.  The  x y,  and  ^r-axis  are  positively  oriented  to  the  east, 
north,  and  “straight-up”  respectively.  The  Earth’s  axis  of  rotation  however,  is  none 
of  these,  but  rather  passes  through  the  North  and  South  Poles.  This  conundrum, 
depicted  in  Figure  2,  is  the  reference  frame  in  which  we  develop  the  geophysical  fluid 
flow  equations. 

Before  continuing  further,  we  revisitthe  two-dimensional  rotating  system  pre¬ 
sented  in  the  previous  section  in  Figure  1.  We  add  a  third  dimension  to  this  system 
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Figure  2.  The  Earth  Model 


by  defining  a  unit  vector  k  at  the  origin  that  is  perpendicular  to  the  plane.  A  rotation 
vector  is  described  as  fl  =  flk.  We  use  this  to  write  the  inertial  velocity  equations 
(11.22)  in  vector  form: 


U  «  —  flu  —flu 

U  =  =  =u+  . 

V  v  -1-  fir  -fflx 


(11.30) 


Now  let  r  =  ri  -f  y'}  and  consider: 

i  J  k 

Jlxr=  0  0  fl  =  fl(-i/i  +  rj).  (11.31) 

x  y  0 

This  allows  us  to  write  the  inertial  velocity  equations  (11.22)  in  condensed  vector 
form: 

U=u+fixf.  (11.32) 

Similarly,  the  inertial  acceleration  equations  (11.28)  in  vector  form  are: 


A  =  a  + 


—2  fir 
2  flu 


(11.33) 
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Using  an  analogous  approach  we  show  that  the  condensed  vector  form  of  the  inertial 
acceleration  equations  is: 


A  =  a  4-  2ft  x  u  4  ft  x  (ft  x  F). 


(H. 34) 


In  this  form,  the  Coriolis  acceleration  is  2ft  x  u  and  the  centrifugal  acceleration  is 
ft  X  (ft  X  r).  With  respect  to  (11.32)  and  (11.34),  the  time  derivative  for  the  inertial 
frame  is  equivalent  to  applying  the  operator: 


D  A 
-)-  ft  x , 
Dt 


(11.35) 


where  —  is  the  derivative  following  a  fluid  introduced  in  (II. 11). 

We  now  consider  two  additional  simplifications  in  the  development  of  the  ac¬ 
celeration  equations  for  the  Earth  model.  First  we  neglect  extraneous  terms  resulting 
from  the  Earth’s  curvature.  In  general  this  can  he  done  if  L  <3C  r  where  r  is  the  ra¬ 
dius  of  the  sphere,  On  Earth,  L  <  1000  kilometers  is  acceptable  [Rjef.  28],  A  second 
simplification  to  (11.34)  comes  from  our  intuition  about  planetary  phenomenon!.  The 
centrifugal  force  is  an  outwardly  normal  force,  or  from  the  vantage  point  of  the  Earth 
observer,  a  force  that  acts  straight  up.  Yet  nothing  ever  is  “flung7’  upward  from  the 
face  of  the  Earth  because  gravitational  forces  keep  centrifugal  forces  in  check.  In  the 
absence  of  rotation,  gravity  would  hold  the  Earth  together  as  a  perfect  sphere.  The 
presence  of  rotation  and  accompanying  centrifugal  forces  distort  the  sphere,  flattening 
it  to  the  extent  that  gravity  and  the  centrifugal  force  negate  each  other.  Hence,  we 
neglect  centrifugal  forces  in  the  Earth  model  and  the  inertial  acceleration  equations 
become: 


A  =  a  —  2ftr, 
B  =  b  -\-  2ftu, 


(11.36) 


or  in  vector  form: 

A  =  a  4-  2ft  x  u. 


(11.37) 


Using  these  equations,  we  continue  with  the  Earth  model. 
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Referring  to  Figure  2  we  write  the  rotation  vector  ft  in  terms  of  the  unit 
vectors  i,  j,  and  k: 

ft  =  ft  cos  4  ft  sdn  ^k,  (11.38) 

where  <fi  is  the  degree  of  latitude  of  the  observer.  Since  i  is  always  orthogonal  to  ft, 
it  does  not  appear  in  the  equation.  Using  (11.35)  we  write: 

Du 


a  =  — —  2 ft  x  u, 

Dt 


where: 


2ft  x  u  = 


i  J  k 

0  2ftcas</>  2ftsin0 


U  V 


w 


or: 


Th  erf  ore: 


Du 

Ox  =  -|-  2ftu>  cos  <f)  —  2ftrsin  <j>, 


(11.39) 


(11.40) 


2ft  X  u  =  (2 ftu;  cos  <j>  —  2ftr  sdn  4>)\  4  2ftu  sin  <ftj  —  2ft«  cos  </>k.  (11.41) 


Dr 

<iy  =  4  2  ftu  sdn  <f>, 


(11.42) 


Dw 

a3  =  — —  —  2ftucosd>. 

Dt 

For  convenience,  we  define  the  following  quantities: 

Coriolis  Parameter:  /  =  2ft  sin 

Reciprocal  Coriolis  Parameter:  fr  =  2ft  cos  0, 
allowing  us  to  rewrite  (11.42)  as: 

Du 

<*x  =  /•»-  fv. 


(11.43) 


— 


Dr  r 

m +  /”' 


(11.44) 


a,  = 


Dw 

~Dt 


-  /„«■ 
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The  angle  <j>,  and  therefore  the  Coriolis  parameter,  is  positive  in  the  Northern  Hemi¬ 
sphere  and  negative  in  the  Southern  Hemisphere.  The  Coriolis  parameter  is  zero  on 
the  Equator  where  <p  =  0.  The  reciprocal  Coriolis  parameter  is  positive  everywhere 
except  at  the  poles  (0  =  ±90)  where  it  is  zero. 

One  final  simplification  is  now  applied.  As  stated  earlier,  our  model  of  geo¬ 
physical  fluid  dynamics  is  based  on  the  shallow  water  assumption  where  the  vertical 
dimension  (D)  is  very  small  relative  to  the  horizontal  dimensions  (£): 

In  other  words,  geophysical  fluid  flow  is  “almost  two- dimensional”  and  therefore  ver¬ 
tical  flow  components  are  negligible  (e.g.  w  <£  u,v  and  a2  ~  0).  Thus  we  rewrite 
(11.44)  as: 

(11.45) 

These  are  the  acceleration  components  of  the  Earth  model  that  are  used  to  generate 
the  momentum  equations  for  the  geophysical  fluid  flow  model. 


Du  , 

Dr  r 

Dt~  fV’ 

av  =  t: — b 

y  Dt~r  J  , 

a3  =  0. 

4.  Forces  Acting  on  a  Fluid  Control  Volume 

Two  types  of  forces  act  on  a  fluid  in  a  control  volume:  Body  forces  which  are 
proportional  to  the  volume  mass  and  surface  forces  which  are  proportional  to  the 
volume  surface  area. 

Gravity,  the  only  applicable  force  in  the  Earth  model,  acts  along  the  z-axis 
toward  the  center  of  the  Earth.  It  acts  on  a  rectangular  volume  element  with  sides 
dx,  dy,  and  dz  and  is  given  by: 


f gravity  =  —mg  =  —  (p  dx  dy  dz)g  =  —pgdV, 


(11.46) 


where  dV  is  the  incremental  volume  of  the  fluid  element. 

Salient  surface  forces  include  pressure  and  viscosity.  Pressure  acts  in  a  direc¬ 
tion  normal  to  the  volume’s  surface  as  shown  in  Figure  3.  The  total  difference  in 
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Figure  3.  Pressure  Effects  on  a  Rectangular  Fluid  Volume 


pressure  along  dx  is: 


VP. 

—  dx. 
ox 


To  obtain  the  total  force  due  to  pressure  in  the  j- direction,  this  expression  is  multi¬ 


plied  by  the  surface  area  (dy  dz)  upon  which  it  acts.  The  result  is: 


^pressure*  =  dz  = 


(11.47) 


This  force  acts  in  the  direction  of  the  negative  gradient  (e.g.  from  an  area  of  high 
pressure  to  an  area  of  low  pressure).  Similarly  the  expressions: 


ET  _  r  JT/ 

pressure-#  — 


r  —  —  —  dV 

■  pressure.  — 


(II  .48) 


describe  the  force  due  to  pressure  in  the  y-  and  z- directions: 


Viscosity  manifests  itself  in  surfaces  forces  called  shear.  Figure  4  summarizes 
the  viscous  stress  tensor.  Viscous  stresses  are  symmetrical  (eg.  =  r^).  A  summa¬ 


tion  of  viscous  forces  in  the  x-  directi  on  yields: 


~^dV  J  dx  dz  + 


dz  dx  dy, 


(11.49) 
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The  surface  forces  and  body  forces  are  combined  to  generate  the  total  force  equations: 


(11.52) 


We  now  have  the  necessary  expressions  to  set  up  the  momentum  equations  for  the 
geophysical  fluid  flow  model. 


5.  The  Momentum  Equation  for  Geophysical  Fluid 
Flow 

Applying  Newton’s  Second  Law  (Ft  =  mat  =  pa^dV)  and  using  (11.45)  and 
(11.52),  we  now  write  the  momentum  equations  for  the  geophysical  flow  model: 


dp  dr^c  drxy  drX3 
dx  dx  ^  dy  ^  dz  ’ 


dp  drgy  ^ryy  dry3 
dy  dx  dy  dz  ’ 


(11.53) 


pa3  — 


dp  drX3  dry3  dr33 

-  -&-«+-a7+iy  +  -a7 


These  equations  are  simplified  using  the  assumptions  described  earlier: 


•  The  Fluid  is  Incompressible :  Density  is  independent  of  pressure,  and  therefore 
the  model  is  uncoupled  from  thermodynamic  considerations. 

•  The  Fluid  is  Inviscid:  All  viscous  forces  are  equal  to  zero. 

•  The  Fluid  zs  Homogeneous :  We  need  not  deal  with  the  complexities  of  density 
stratification  (this  simplification  will  be  relaxed  in  Chapter  V). 

•  Centrifugal  Forces  are  Negated  by  Gravity.  This  allows  simplification  of  the 
acceleration  terms. 
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•  L  r:  This  allows  us  to  neglect  teims  that  arise  in  the  acceleration  compo¬ 
nents  that  result  from  the  curvature  of  the  Earth  (e.g.  see  Cushman- Rjoisin 

(1994)  [Ref.  28]). 

•  The  Flow  w  Primarily  Horizontal :  Velocity  and  acceleration  terms  in  the  z- 
directdon  are  negligible. 


Since  viscosity  terms  are  negligible  (11.53)  becomes: 


(11.54) 


n  dP 
a=---pg. 


D 


Expanding  the  operator  —  introduced  in  (II.  11)  produces  the  momentum  equations 
in  simplified  form: 


i-momentum: 


du  ,  du  ,  du  .  1  dp 


y-  momentum: 


dv  dv  dv  1  dp 

oT  +  «  oT  +  r  oT  +  /«  = 
at  ox  ay  p  ay 


(11.55) 


^-momentum: 


n  13 'P 

0  =  —r--g. 

pdz 


Since  geophysical  flow  is  “primarily  horizontal7’,  the  term  w—  does  not  appear  in 

az 

the  x-  and  ^-momentum  equations.  Equation  (11.55)  together  with  the  continuity 
equation  (II. 6)  are  used  to  construct  a  shallow  water  model. 


6.  Governing  Equations  for  the  Shallow  Water  Model 

The  shallow  water  model  is  depicted  in  Figure  5.  We  define  a  variable  h  as 
the  height  of  the  fluid  above  a  reference  level  z  =  0.  The  variable  h  varies  with  x, 
y,  and  t  and  describes  the  fluid  action  at  the  surface.  The  “bottom”  of  the  shallow 
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Figure  5.  The  Shallow  Water  Model 


water  environment  consists  of  a  rigid  suiface  whose  height  that  does  not  vary  with 
time  and  is  given  by  z  =  hB(x,y).  We  let  H  =  h(x,  y,t)  —  hB(x,y)  where  H  is  the 
depth  of  the  fluid  layer  with  respect  to  the  bottom  contour.  If  L  is  the  horizontal 
scale,  then  by  the  shallow  water  assumption  it  is  true  that  H  C  L. 

It  is  now  possible  manipulate  the  governing  equations  for  the  Earth  model  to 
produce  the  shallow  water  model.  Integrating  the  z-component  of  the  momentum 
equation  (11.55)  yields: 


p(x,y,z,t)  =  —pgz  -I-  p(x,  y,  t)  (II .56) 

On  the  surface,  z  =  h(x,y,t),  pressure  p  equals  some  ambient  pressure  P0.  Therefore: 

p(x,  y,  z,  t)  =  pg[k(x,y,t)~  z]  +  P0.  (11.57) 

Dropping  the  variable  dependencies  for  brevily,  (11.57)  is  rewritten  as: 

p  =  pg(k  -  z) -\- P0.  (11.58) 


From  (11.58)  we  obtain  pressure  gradients  in  the  horizontal  x  and  y  directions: 


dp 

dx 


dp  dh 

9y=P%' 


(11.59) 
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Thus  we  rewrite  the  x-  and  y-momentum.  expressions  (11.55)  as: 


x-  momentum: 


y-momentum: 


8u  du  du  dh 

T7  +  “5;  +  V  AT,  -  fV  =  ~g^r' 
Ot  Ox  Oy  Ox 


dr 

dt 


dr  dr 

+  «  +  l' T  /«  = 

Ox  Oy 


(11.60) 


Equation  11.59  also  implies  that  the  horizontal  acceleration  components,  a x  and 
dy,  given  by  (11.45)  and  (11.55)  are  independent  of  z.  It  follows  that  the  horizontal 
velocity  components,  v.  and  t,  are  independent  of  z.  This  allows  us  to  uncouple  (II. 6) 
and  solve  for  w(x,y,  z,t ): 


w{x,y,z,t)  =  -z 


dx 


dy 


+  w(x,y,t) 


(11.61) 


Now  consider  the  flow  along  the  bottom  contour  hB{x,  y).  Since  the  contour  is  rigid, 
there  is  no  normal  flow,  Any  velocity  in  the  z-  direction  is  due  to  fluid  flowing  tan¬ 
gent  to  the  bottom  contour.  Figure  6  depicts  this  consideration  for  the  x—  and  z— 
directions.  If  u  is  the  velocity  component  in  the  x- directi  on  flowing  tangent  to  the 


Figure  6.  Flow  Along  the  Bottom  Contour  in  the  j- Direction 


bottom  contour,  its  contribution  to  in  is: 


v.  tan(a) 


u- 


dhB 

dx  ’ 


(11.62) 
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where  a  is  the  angle  of  the  line  tangent  to  hB  and  x-axis.  Similarly  in  the  y- direction, 


the  contribution  to  w  is; 


v  tan  (if)  =  r 


BhB 

~dy’ 


(11,63) 


where  0  is  the  angle  of  the  line  tangent  to  hB  and  y-axis.  It  follows  then  that: 


,  t  dkB  ,  dhB 

w{x,y,z  =  hB,t)  = 


(11,64) 


Therefore  from.  (11,61)  and  (11,64)  we  generate  an  expression  for  w  along  the  bottom 
contour : 

{ fin  fii'  \  fih  fih  t> 

(11,65) 


.  /  Bu  <9r\  dhB  dhB 
v)  =  {hB  -  z)  —  +  3-  +  «—  4- 1- 


KBx  By/  Bx  By  ' 

A  corresponding  condition  at  the  surface,  z  =  h  yields: 

Bh  Bh  Bh 

w{x,y,z=  h,t)  = —  -\-u--\-v-.  (11,66) 

Bh  .  .  , 

Since  the  surface  is  not  rigid,  we  pick  up  in  the  expression.  Therefore,  (11.65)  and 

at 

(11.66)  generate: 


Bh  Bh  Bh 

a+”al  +  t^ 


’  Bu  Bv 
kBx  +  By 


Bhi 


Bhi 


Bx 


By’ 


(11,67) 


which  is  simplified  to: 

(11.63) 

Equations  (11.60)  and  (11.63)  are  the  governing  equations  for  the  shallow  water  model. 
In  its  current  form  the  model  is  non-linear.  We  desire  a  linear  form  for  subsequent 
investigations. 


dh  d  d 

Bt+B^^l('k~  +  By  ^  “  ^b)]  =  °‘ 


7.  Linearizing  the  Shallow  Water  Model 

To  linearize  the  shallow  water  model  we  conduct  a  perturbation  analysis  on  the 
x-  and  y- momentum  equations  (11.60)  and  the  vertical  momentum  equation  (11.63). 
Perturbation  analysis  on  the  horizontal  flow  is  accomplished  by  assuming  that  the 
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Figure  7.  Shallow  Water  Model  Cross  Section  with  h  =  H  -f  hs 


x-  and  y- components  of  velocity  are  dominated  by  constant  terms  ( U  and  V’).  Su¬ 
perimposed  on  these  are  small  variations,  u*(x,y,  t)  and  r*'(x,  y,  t).  Mathematically 
stated: 

u  =  U  -f  u*,  where  V  is  constant  and  ir*  ~  0(6)  <£  U, 
t  =  V  -fr*,  where  V  is  constant  and  r*  ~  0(6)  <£  V. 

Applying  these  to  (11.60)  yields: 


du*  du*  du *  dk 

_  +  (t;+K.)_+(V-+r.)__/(v-  +  ,.)  = 


Si'*  .  -T-  ..Si'*  r-TT 

ST  +  ^+^)-^  +  {V+v)-^  +  f{V  +  v)  =  -g-. 
Ignoring  terms  of  0(52)  results  in  the  following  simplification: 


du*  Tjdu*  _du*  _¥_  dk 

+  - /(V  +  r  )  =  -g 

dy 


dt 


dx 


dx  ’ 


dv*  TJdr*  T.dv*  ,  dh 

ar+t,aT  + v  n; +  w  + u )  =  -%■ 


(11.69) 


(11.70) 


(11.71) 


A  similar  perturbation  approach  is  applied  to  the  vertical  momentum  equation 
of  the  shallow- water  model  (11.68).  We  let  k( x,  y,  t)  =  kB(x,  y)-\-H(x,  y,  t)  as  depicted 
in  Figure  7,  Applying  this  to  (11.68)  yields: 
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where  H  =  h  —  kg.  We  refine  H(x,  y,  t )  further  into  two  parts  as  depicted  in  Figure  8. 
The  component,  Hq(x,  y)  is  resultant  from  the  constant  velocity  terms  V  and  V 


Figure  8,  Shallow  Water  Model  Cross  Section  with  h  =  H0  17  -f  hB 

introduced  in  (11.69).  If  U  and  V  are  zero  then  H0  is  constant.  Likewise,  if  the 
bottom  is  fiat,  then  there  is  no  velocity  in  the  vertical  direction  along  the  bottom 
that  results  from  a  non-zero  U  and  V,  and  again  H0  is  constant.  Superimposed  on 
Hq  is  a  small  amplitudinal  variation  that  represents  the  wave  action  and  is  given  by 
rj(x,y,t).  Mathematically  stated  this  is: 

H(x,y,t)  =  H0(x,y )  -+  rj(x,y,t)  where  77  ~  0(5)  <  H0  (11.73) 

Applying  these  terms,  as  well  as  the  perturbation  terms  introduced  for  the  horizontal 
velocity  allows  us  to  rewrite  (11.72)  as  follows: 

f[  +  £  [(£>  +  «*)(»«  +  17)]  +  ^  [O'  +  **)(«>  +  >7)1  =  0-  (11.74) 

Neglecting  terms  of  0(52)  gives  us: 

^  \U(H0  + ,)  +  «•  H0]  +  {V(H0  +  ,)  +  V‘H„]  =  0. 


(11.75) 


(11.76) 


Since  U  and  V  are  constant,  11.75  becomes: 

o  d  d  d  d 

~dt  ^  ~dx  ^  h  ^  £^(#o  4-  -f  ^^(t^i/o)  = 

This  can  be  further  simplified  if  the  bottom  is  flat  (e.g,  H0  is  constant): 


(11.77) 


Finally,  if  we  assume  that  V  and  U  are  zero  (e.g.  no  advectdon),  the  equation  becomes: 

(II. 

This  is  the  simplest  form  of  the  linearized  vertical  motion  component  of  the  shallow 
water  model. 

We  now  revisit  the  horizontal  flow  equations  (11.71).  By  including  the  refine¬ 
ments  and  perturbations  on  h,  (11.71)  becomes: 


<9u*  TIdu1'  T„du*  ( dfiB  dH0  dn 

i>r+t,ar+vftr-*v+l'>  =  -Har  +  ^r+Sr 


%  +  v%  +  v%- -  -» 


dhB  dHo  drt ] 

~di  +  ~di  +  d^ 


(11.79) 


Since  we  have  assumed  no  advectdon  (Z7,  V  =  0  and  therefore  Hq  is  constant)  and  that 
the  bottom  is  flat  ( hB  =  0),  the  linearized  horizontal  flow  component  of  the  shallow 
water  equation  can  be  stated  as: 

(11.80) 

Equations  11.78  and  11.80  are  the  governing  equations  for  the  linearized  shallow  water 
model.  They  will  be  further  refined  in  the  following  section. 
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8.  Deriving  the  Klein-Gordon  Equation 

We  now  combine  linearized  components  of  the  shallow  water  model  to  obtain 
the  well  known  Klein- Gordon,  oi  dispersive  wave  equation.  A  step-by-step  derivation 
follows. 

Step  1:  Perform  the  following  substitutions  to  (11.78)  and  (11.80): 


u  =  v*Hq  and  i  =  v  *H0. 


The  result  is: 


-  -  -  -gH*-, 


&  ,  rr&l 

_  +  /„  =  -gHa- 

dndidi 

ai+ai  +  a^  =  0- 


(11.81) 

(11.82) 

(11.83) 

Step  2:  Assume  /  is  constant  and  take  the  partial  derivative  of  (11.81)  with  respect 
to  r: 

A  (^\  _f*  =  *2 

dx  J  dx  °  dx 2’ 

and  the  partial  derivative  of  (11.82)  with  respect  to  y : 

d  (Bv\  .di  _  d2v 

dy{dt)+fdy~  gH°dy>’ 

and  add  the  resulting  equations: 

d  f  di  <9t\  /  <9r  <96  \  2 

+  =  ~gH° V  V ' 

Step  3:  Take  the  partial  derivative  of  (11.81)  with  respect  to  y\ 

d  /  di\  di  <9  /  dxj  \ 

dy  5  °dy  V<9xJ  ’ 

and  the  partial  derivative  of  (11.82)  with  respect  to  x: 

d  ( di\  di  d  ( dn\ 

di  ( at )  +  fdi  =  ~gH°di  ( dU)  ’ 


(11.84) 
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ajid  subtract  the  resulting  equations: 

d  ( di  di 


d  /  di  di  \  /  <96  \ 

<9i  l  <9*/  <9x)  t  <9y  <9x  J 


(11,85) 


Step  4:  Take  the  partial  derivative  of  (11,84)  with  respect  to  t  and  rearrange  the 


resulting  terms: 


„  d  f  dv  di  i  cr  1  Ou  crv  \  „  d  \ 

f  dt  (  d^  ~  =  m  (  d^  +  d^ )  "*■  gH°dt  (V  7?) ' 


d 2  (di  di 


|_<9t  \dx  dy  j  J  dt2  \<9x  dy  j  dt 

Multiply  (11.85)  by  /  and  rearrange  the  resulting  terms: 


d  f  di  di 
dt  l  <9x  dy 


di  |  <96  \ 

dy  dx)  ' 


Summing  these  yields: 


di  di 


[ae  +  fj  +  dy)  +  gHodt  O7^)  ~  °- 

Step  5:  Using  (11,83)  rewrite  (11,86)  as: 

Step  0:  Integrating  (11.87)  with  respect  to  t  and  letting  Cq  =  gHo  yields: 


(11,86) 


(11,87) 


^-C^V+fV=S(x,y), 


(11,88) 


where  S(x,  y)  is  an  arbritrarv  function  of  intergration.  Equation  11,88,  the  linear 
inhomogeneous  form  of  the  Klein-Gordon  equation,  is  a  restatement  of  the  linearized 
shallow  water  equation,  It  also  describes  other  behaviors  such  as  lateral  vibrations 
of  membrane  strips  and  acoustic  pressure  waves  in  dispersive  media  [Ref,  24],  We 
continue  with  analytical  considerations  of  the  Klein-Gordon  equation  and  the  concept 
of  dispersive  waves  in  the  next  section. 

9.  Analytic  Considerations  of  the  Klein-Gordon 
Equation 

Consider  a  homogeneous  form  of  the  Klein-Gordon  Equation  in  one  dimension: 


&V  _  C2  A  ,  f2  =  0 

dt2  c° dx2+fv 


(11,89) 
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This  equation  has  the  following  solution: 


=  exp[z(fcx  —  ut)], 


(11.90) 


where: 

^2=C2fc2-f/2.  (II.91) 

Here  to  is  the  angular  frequency  and  k  is  the  wave  number,  Solving  for  k  one  obtains 
the  dispersion  equation  for  the  Klein- Gordon  equation: 

vtM5 


k  = 


(11.92) 


The  phase  velocity  of  the  wave  (e.g.  the  speed  of  a  wave  crest)  is  given  by: 


to 

tp  =  T 


(11.93) 


Group  velocity  is  the  velocity  at  which  the  wave  energy  propagates  and  is  given  by: 

dto 


l'9  dk 


(11.94) 


If  the  phase  and  group  velocities  are  not  equal,  the  wave  is  dispersive  and  the  wave 
shape  deforms  as  it  travels.  Using  (11.91)  we  have: 


\J Gq  fc2  T  /2 
k 


= 


kCl 


jQV  +  p 


(11.95) 


If  /  ^  0,  then  ip  ^  rfl,  and  the  solution  to  the  Klein- Gordon  equation  produces 
dispersive  waves.  If  /  =  0,  then  vp  =  rfl  =  Co  and  the  resulting  wave  is  non- 
dispersive.  With  respect  to  the  Earth  model,  /  is  the  Coriolis  parameter  given  by 
(11.43).  The  magnitude  of  /  increases  as  we  go  north  or  south  from  the  equator,  Thus 
dispersion  effects  will  increase  away  from  the  equator.  Hence  the  rotating  Earth,  with 
the  exception  of  the  equator  where  /  =  0,  is  a  dispersive  environment  with  respect 
to  the  shallow- water  model. 

This  form  of  the  Klein- Gordon  equation  is  used  in  initial  investigations  of  the 
Higdon  NRBC  because  it  is  a  relatively  simple  mathematical  model  of  dispersive  wave 
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behavior.  These  effects  are  an  integral  partin  the  development  of  Higdon  NRBC’s. 
The  simplicity  of  the  equation  also  makes  possible  comparisons  between  numerical  and 
analytical  solutions  and  provides  for  easy  testing  of  proposed  boundary  conditions. 
Chapter  2  follows  with  a  description  of  the  Higdon  NRBC. 
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III.  HIGDON  NON-REFLECTING 
BOUNDARY  CONDITIONS 

Our  goal  is  to  accurately  describe  the  propagation  of  dispersive  waves  in  a 
“manageable7’  domain.  Such  a  domain  is  one  that  will  not  overwhelm  computer 
capabilities  when  applying  numerical  techniques  to  its  interior  and  boundaries.  The 
actual  domain  in  which  the  wave  travels  is  much  larger,  or  perhaps  infinite.  To  restrict 
ourselves  to  a  smaller  domain  of  interest,  artificial  boundaries  are  constructed  that 
allow  waves  that  impinge  upon  them  to  freely  pass. 

Constructing  the  mathematical  analog  for  a  NRBC,  is  elusive.  Schemes  exist 
that  allow  for  the  total  absorption  of  non- dispersive  waves.  However,  when  dealing 
with  dispersive  waves,  the  results  are  less  than  perfect,  and  spurious  reflections  occur 
at  the  artificial  boundary.  Hence  the  NRBC  problem  is  one  of  optimization  (e.g.  min- 
imiang  unwanted  reflection  thus  allowing  most  of  the  wave’s  energy  to  pass).  The 
Higdon  NRBC  is  the  focus  for  our  investigation  because  it  exhibits  numerous  advan¬ 
tages.  We  initially  use  a  model  that  requires  a  single  artificial  boundary.  From  this 
we  develop  the  details  of  the  Higdon  NRBC. 

A.  THE  SEMI-INFINITE  WAVE  GUIDE  PROBLEM 

The  semi- infinite  wave  guide  problem  provides  a  vehicle  to  explore  the  proper¬ 
ties  of  the  Higdon  NRBC  and  is  depicted  in  Figure  9.  A  Cartesian  coordinate  system 


r, 

Figure  9.  Semi- Infinite  Wave  Guide 
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(x,  y )  is  introduced  such  that  the  wave-guide  is  parallel  to  the  x- direction.  The  width 
of  the  wave-guide  is  denoted  by  b.  In  the  wave-guide  we  consider  the  inhomogeneous 
Klein-  Gordon  equation  (11.88): 

§£-csw+a-s 

where,  rf  is  the  unknown  wave  held,  Q,  is  the  given  reference  wave  speed,  /  isthe  given 
dispersion  parameter,  and  S  is  a  given  wave  source  function.  C0  and  /  are  functions 
of  location,  but  it  is  assumed  that  outside  a  finite  region  they  do  not  depend  on  x. 
The  wave  source  5  is  a  function  of  location  and  time,  but  is  assumed  to  have  local 
support. 

Dirichlet  or  Neumann  boundary  conditions  are  specified  on  the  north  and 
south  boundaries,  and  17*?: 

>7=0  o,  |=0.  (III.1) 

In  acoustics,  these  correspond  to  “soft  walT  and  ‘hard  wall”  conditions,  respectively. 
A  Dirichlet  boundary  condition  is  prescribed  on  the  west  boundary  T K. : 


>7(0  ,!/,*)  =  %-{y,*),  (in. 2) 

where  rjw(y,t)  is  an  incoming  wave  function.  As  x  — *  oo,  the  solution  is  bounded 
and  does  not  include  incoming  waves.  To  complete  the  problem  statement,  the  initial 
conditions: 

fhi 

y,  0)  =  rjo,  ^(x,y,0)  =  wo,  (III.3) 

are  given  at  time  t  =  0  for  the  entire  domain,  We  assume  that  the  functions  7j0  and 
wo  have  local  support. 

We  now  truncate  the  semi-infinite  domain  by  introducing  an  artificial  east 
boundary  B  at  x  =  xr  which  we  call  De.  This  boundary  divides  the  original  semi- 
infinite  domain  into  two  subdomains:  an  exterior  domain  V,  and  a  finite  computa¬ 
tional  domain  bounded  by  Tjv,  Ts,  De,  and  T h-  .  We  chose  the  location  of 
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such  that  the  entire  support  of  5,  rj o,  and  wq  is  in  ft.  Thus  on  Tf  the  homogeneous 
counterpart  of  the  Klein- Gordon  equation  holds: 

§j?-<Sv«,+A- 

On  r ^7,  Cq  and  / 2  are  y- dependent  (or,  as  a  special  case  constant). 


B.  THE  HIGDON  NON-REFLECTING  BOUNDARY 
CONDITION 

The  Higdon  condition  for  non- dispersive  acoustic  and  elastic  waves  was  pre¬ 
sented  and  analyzed  in  a  sequence  of  papers  (see  e.g.  [Rjef.  17]  -  [Ref.  20]).  These  were 
later  extended  to  include  dispersive  waves  [Ref.  21].  To  obtain  a  well- posed  problem 
for  the  finite  domain  ft  (see  Figure  9)  we  impose  a  reformulation  of  the  Higdon  NRBC 
[Ref.  21]  on  Re  to  reduce  spurious  wave  refection. 

The  Higdon  NRBC  is  obtained  by  composing  simple  first-order  differential 
operators.  For  example,  a  Higdon  NRBC  of  order  J  (denoted  by  Hj)  is: 


(III  .4) 


In  this  section  we  will  show  that  Higdon  NRBC’s  have  several  ad  vantages  including; 

•  Their  reflection  coefficients  are  easily  determined. 

•  They  are  exact  for  all  waves  that  propagate  in  an  x-direction  with  phase  speeds 
equal  to  either  of  Ci  through  Cj. 

•  They  constitute  a  sequence  of  conditions  of  increasing  order.  No  asymptotic 
approximation  is  involved  in  their  construction,  enabling  one  to  obtain  solu¬ 
tions  with  unlimited  accuracy. 

•  They  are  robust.  Reflection  coefficients  become  smaller  as  the  order  J  in¬ 
creases.  A  good  choice  of  C)’s  leads  to  better  accuracy  for  a  smaller  J,  but 
reductions  in  spurious  reflection  can  still  be  obtained  with  non-optimal  C/s 
by  simply  increasing  J . 

•  They  are  very  general  applying  to  a  variety  of  wave  problems  in  one,  two, 
or  three-dimensional  configurations.  Moreover,  they  can  be  used  for  wave 
problems  in  dispersive  and  stratified  media. 
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To  discover  the  motivation  for  the  form  introduced  above  (III .4),  consider  a 
possible  solution  to  the  semi-infinite  wave  guide  program  suggested  by  (11.90): 


rj  =  AYn{y)  exp[zfc(x  -  Cxt  +  <■)], 


(III. 5) 


where  A  is  the  wave  amplitude,  %•  is  the  phase,  k  is  the  x- component  wave  number, 


u? 


u>  is  the  wave  frequency,  and  Cx  =  —  is  the  phase  velocity.  Yn{y)  is  determined  from 

nr 

the  dependency  of  Co  and  /  on  y.  Rjecall  from  the  Klein  Gordon  equation  derivation 
that  Co  =  gH o,  where  Ho(x,y)  is  described  in  Figure  8.  Since  this  analysis  assumes 
a  flat  bottom,  H0,  and  therefore  Co,  is  constant.  Also  recall  that  /,  the  dispersion 
or  Coriolis  parameter,  varied  with  the  latitude  (11.43).  With  respect  to  the  Earth 
model  (Figure  2),  <j>  was  determined  to  be  a  function  of  y  only.  Therefore,  in  the 
shallow  water  model,  /  is  a  function  of  y  only.  Now  consider  the  real  part  of  (III.  5): 


rj  =  AYn{y)  cos[fc(x  -  Cxt  -f  if')]- 


(III. 6) 


Substitution  into  (HI .4)  generates: 


AV^(i/)fcsin  [fc(x  —  Cxt  V')] 


n  (c.  -  c,) 


=  0, 


(III  .7) 


which  implies: 

n(C*  -  Cj)  =  0.  (HI .8) 

Thus  the  Higdon  NRBC  is  exact  (eg.  no  portion  of  the  wave  is  reflected)  at  He  if  the 
phase  speed  Cx  matches  any  of  the  chosen  Higdon  parameters  Cj . 

Typically  solutions  to  the  Klein- Gordon  equation  consist  of  an  infinite  number 
of  waves,  whereas,  we  are  limited  to  the  selection  of  a  finite  number  of  parameters. 
Therefore,  in  most  cases  Cx  ^  Cj.  We  can  still  validate  (III. 8)  by  assuming  that  the 
impinging  wave  splits  at  rr  into  a  reflected  and  passing  wave.  The  magnitude  of 
the  reflected  wave  is  easily  analyzed.  Consider  a  simplified  form  of  (III. 6)  in  which 
Yn{y)  =  1  and  =  0: 

rj  =  Acos[£(x  —  Cxt)],  (III. 9) 
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ajid  a  first-order  Higdon  NRBC: 


(111,10) 


where  Cx  ^  Cj.  The  original  wave  impinging  on  Tr  is  the  sum  of  the  wave  that 
passes  through  and  the  wave  that  reflects  hack  into  the  domain,  Mathematic  ally 
stated: 


A  cos[fc(x  —  C^f)]  =  Ar  cos[— k(z  -1-  Cxi)\  4-  Ap  cos[fc( j  —  Cai)].  (Ill, 11) 

where  Ar  and  Ap  are  the  amplitudes  of  the  reflected  and  passing  waves  respectively. 
Note  that  reflection  affects  the  wave  hv  reversing  its  direction  of  travel,  or  mathe¬ 
matically  by  reversing  the  sign  of  the  wave  number  k.  The  wave  frequency  to  for 
the  reflected  wave  remains  unchanged.  If  any  reflection  occurs,  then  0  <  p4r|  <  \A\. 
With  regards  to  the  passing  wave,  the  wave  number  k  and  phase  speed  Cx  remain 
unchanged.  However,  \Ap\  will  be  reduced  with  respect  to  \A\.  Substituting  the 
right-hand  side  of  (III .11)  into  (III.  10)  yields: 

Ar(Cx  A  Cj )  -  AP(CX  -  C})  =  0.  (III. 12) 


Using  this  equation  we  define  the  reflection  coefficient  R  for  a  first-order  Higdon 


NRBC  to  be: 


Ar 

Cj-Cx 

Ap 

Cj-f  c* 

(111,13) 


This  yields: 

KI  =  ITflMI  “d  ia,i  =  TTbw  <III14> 

We  see  from  (III.  13)  and  (III  .14)  that  |A-|  <  |-4p|  when  Cx  ^  0.  Note  that  R  — ?•  1  as 


Cjx  — ^  0  implying  that  waves  with  low  phase  speeds  will  result  in  maximum  reflection 
of  |^lr|  =  \AP\  =  .5 |j4|  at  the  artificial  boundary.  This  circumstance  is  mitigated  in 
that  these  same  low-speed  waves  reach  the  boundary  at  times  that  are  often  outside 
the  scope  of  the  problem. 
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The  reflection  coefficient  for  a  J^-order  Higdon  NRBC  is  easily  derived  using 
the  above  techniques: 

j 

«- n 

j=i 


Cj-Ca 


Cj  +  Cx 


(III. 15) 


This  expression  represents  the  product  of  J  factors  that  are  less  than  one.  Hence  sim¬ 
ply  increasing  the  order  J  reduces  the  amplitude  of  the  reflecting  wave.  Theoretically, 
one  could  reduce  the  amplitude  of  the  reflected  wave  at  r£  to  zero  without  regard  for 
the  value  of  Cj  by  letting  J  —*■  co.  Practically,  however,  we  are  limited  by  computer 
capability.  We  can  only  select  a  finite  number  of  Cj’s  and  therefore  must  tolerate 
some  spurious  reflection.  Fortunately,  R  can  by  reduced  significantly  by  intelligently 
selecting  Cj’s.  Strategies  for  this  are  discussed  in  the  next  section. 


C.  DETERMINING  OPTIMAL  VALUES  FOR  THE 
HIGDON  PARAMETERS 

Existing  literature  offers  no  analytical  means  to  optimize  the  value  of  the  Hig¬ 
don  parameters  for  a  dispersive  wave.  Three  general  methods  have  been  suggested: 
(1)  a-priori  selection,  (2)  computer  automated  selection,  and  (3)  dynamic  selection. 
The  first  two  methods  are  explored  in  this  section. 

The  first  general  method  was  suggested  by  Higdon  (see  e.g.  [Rjef.  17]  -  [Ref. 
21])  and  selects  Cj’s  using  an  “educated  guess75 .  One  takes  advantage  of  the  properties 
of  the  reflection  coefficient  R  and  parameters  of  the  Klein- Gordon  equation  to  accom¬ 
plish  this.  A  second  general  method  chooses  Cj’s  automatically  by  computer  code 
as  a  preprocess.  These  methods  typically  use  information  about  the  interior  wave  to 
select  Cj’s  that  minimize  reflection.  Givoli  and  Neta  suggest  a  simple  approach  in 
which  Cj’s  are  determined  from  wave  numbers  that  are  evenly  distributed  over  the 
span  of  maximum  resolvable  wave  numbers.  In  another  approach  they  recommend 
using  the  minimax  method  to  pick  k  values  [Ref.  24] .  To  test  these  methods,  the 
“oscilloscope75  method  is  developed  to  fine  tune  the  Cj’s  and  minimize  the  reflection 
of  a  known  wave.  Theoretically,  this  procedure  produces  the  best  result,  but  it  is  too 
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expensive  computationally  to  seive  as  a  preprocess.  It  is,  however,  a  useful  measure 
of  the  effectiveness  of  other  suggested  methods. 

To  quantitatively  compare  the  different  schemes,  several  numerical  experi¬ 
ments  are  conducted  for  a  semi-infinite  one- dimensional  domain  [0,  co)  on  the  x-axis. 
This  domain  allows  waves  to  travel  in  the  x-directdon  only,  eliminating  a  need  to  con¬ 
sider  geometric  dispersion  which  occurs  in  the  semi-infinite  channel.  It  also  allows  us 
to  determine  analytically  the  effects  of  Higdon  parameters  Cj  on  the  reflection  coeffi¬ 
cient  H.  The  waves  in  the  domain  are  governed  by  the  one- dimensional  homogeneous 
Klein- Gordon  equation: 

!?-<4?+a=o.  (hi. i6) 

An  artificial  boundary  is  placed  at  4 tt.  A  continuous  wave  is  assumed  to  exist 
inside  the  truncated  domain  with  the  following  initial  conditions: 

76  1 

rj(z,  0)  =  —  cos(fcx)  and  rjt(z,  0)  =  0.  (III. 17) 

fc=L  k 

The  solution  is  assumed  to  be  of  the  form: 

75 

rj(z,t)  =  ^  Afccos(fcx  —  vt),  (III. 18) 

te=L 

where  the  dispersion  relation: 

u;2  =  C02fc2  +  /2  (III. 19) 

is  necessary  to  satisfy  (III.  16),  Substituting  (III. 19)  into  (III. 18)  yields: 

J7(x,  f)  =  Yi  ^2  cos  (fcr  “  \j Cok2  T  f2  ■  (III. 20) 

The  resulting  dispersive  wave  is  somewhat  contrived  since  the  number  of  wave  num¬ 
bers  k  is  finite.  However,  it  serves  for  the  current  purpose  to  analyze  various  schemes 
to  select  Higdon  parameters,  Rewriting  (III. 20)  generates: 

*)  =  Y2  h  cos[fc(' 1  ~  C k*)L  (In  '21) 

fc=L  K 
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where: 


Cfc 


^/cgfc2  4-  P 
k 


(III. 22) 


Note  that  C*  >  Co,  a  fact  that  is  important  later  when  estimating  Cj’s.  Equation 
(III  .21)  is  used  to  determine  the  reflection  coefficient  R : 


R  = 


Cj-Cfc 

Cj  T  Cfc 


J 


=  n 


fcCj  -  ^jciv  +  p 
kCi+yjClkl+P 


(III. 23) 


Thus  the  equation  for  the  reflected  wave  is: 

t  P  (flfr)  c«<-fcr  -  Jose  +  f‘  t ),  (III .24) 

which  can  "be  use  to  determine  the  error  || e(t)  ||  resulting  from  a  specified  Higdon 
NRBC.  Note  that,  by  (III. 23),  R  is  a  function  of  k  and  must  appear  inside  the 
summation.  We  calculate  the  error  by  taking  the  norm  of  the  reflected  wave  from  0 
to  Tr.  Since  the  experimental  data  is  discrete,  the  equation  for  the  2- norm  of  the 
reflected  wave  at  time  t  is: 


won = e  (hi  .25) 

*=i  ' 

where  Nx,  the  number  of  elements  in  the  x- vector,  is  determined  by  the  refinement 
of  the  grid.  In  each  experiment  we  set  Co  =  2  and  /  =  1.  The  error  measurement  is 
taken  at  time  t  =  100  time  units. 


1-  Experiment  One:  =  C0 

In  the  first  experiment,  we  disregard  the  wave  generated  inside  the  domain 
and  offer  our  best  guess  for  determining  Cj’s.  As  shown  by  (III. 22),  Cfc  >  C0.  Fur¬ 
thermore,  from  (III. 23),  we  see  that  for  |C0fc|  »  |/|, 

(III. 26) 

It  is  therefore  logical  to  let  Cj  =  Co  for  j  =  1...  J.  For  a  first  order  Higdon 
NRBC  (Hi),  this  is  equivalent  to  the  Sommerfeld  condition  which  is  utilized  widely 


0=1 


Cj-Cp 

Cj  TCo 
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in  meteorology  [FLef .  24],  By  simply  increasing  the  order  J  of  the  Higdon  NRBC, 
(III  .15)  revealed  that  the  reflection  coefficient  R  is  reduced.  When  we  let  C0= 2  and 
/=  1,  our  example  wave  (III. 20)  becomes: 


75 


r/(x>  *)  =  51  Tj  cos(fcx  -  %/4 k2  -f  1  t). 

fc=L  k 

The  Jth- order  reflection  coefficient  (III. 23)  for  this  wave  is: 


(III. 27) 


m  =  n 


2k  -  V4 k2  +  1 


(III. 28) 


2k  +  V4fc2  +  1 

Using  this  equation  we  generate  the  reflected  wave  resulting  from  the  boundary  con¬ 
dition.  Figure  10  (left)  displays  the  reflected  wave  for  a  first-order  (i/L)  through 
third-order  (H2)  Higdon  NRBC.  For  H 1  through  Hg  the  amplitude  of  the  reflected 
wave  has  decreased  to  the  point  that  it  is  not  visible  unless  a  smaller  scale  is  used. 
Figure  10  (right)  displays  the  norms  of  the  reflected  waves,  which  are  taken  to  be  the 
measurement  of  error  for  the  boundary  condition.  This  figure  quantifies  the  rapid 
drop  off  in  the  reflected  or  error  norm  for  increasing  J. 


Figure  10,  Left:  Experiment  la,  at  t =100 :  Cj  =  C0  with  C0=2  and  f=l  (Solid 
Line  Depicts  Hi.  Dotted  Line  Depicts  H2  Plot).  Right:  ||  at  t=100  for  Hi  through 

Hq. 


We  might  refine  this  method  even  further  by  considering  the  following  argu¬ 
ment:  Although  (III.  23)  reveals  Cj  — *  C0  as  k  gets  large,  (III. 28)  suggests  that  R  is 
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small  for  large  k  anyway.  Therefore  "by  selecting  Cj  =  Cb,  we  select  Higdon  coeffi¬ 
cients  that  damp  out  waves  that  had  a  low  lelection  coefficient  R,  while  waves  with 
low  wave  numbers,  and  therefore  higher  reflection  coefficients,  are  not  considered  in 
the  estimation  of  Cj.  Therefore,  using  (III. 22),  we  consider: 


(III. 29) 


Improved  results  are  reported  in  Figure  11.  However,  one  has  presupposed  some 


knowledge  of  the  wave  action  inside  the  finite  domain  in  order  to  estimate 


2.  Experiment  Two:  Cj  Determined  from  Wave  Num¬ 
bers  Distributed  Evenly  over  [Amin,  kmax] 

In  the  second  experiment,  we  selected  C/s  "by  evenly  distributing  wave  num¬ 
bers  over  the  interval  k^a^]  ■  All  other  parameters  used  are  the  same  as  those 

used  in  the  first  experiment.  Here,  fc^tn,  the  minimum  significant  wave  number,  might 
be  determined  by  some  advanced  knowledge  of  the  internal  wave,  or  by  the  time  scale 
of  the  problem  (e.g.  k^^  describes  a  wave  that  moves  so  slowly  that  it  never  reaches 
the  boundary).  For  this  experiment,  k,,,^  =  1.  With  regards  to  k max,  assuming  10 
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grid  points  pei  wave  length  are  necessary  to  resolve  any  given  wave,  a  reasonable 
estimate  is: 

fcf m*  =  (III. 30) 

where  Aj  is  the  grid  spacing.  The  k  values  are  then  given  by: 

kJs  =  i  =  0, 1, ...  J  -  1  (HI  .31) 


where  J  is  the  order  of  the  Higdon  NRBC  and  kj.  is  the  3™  fc-value  for  a  J™-order 
boundary  condition.  Equation  (III. 30)  suggests  that  kma*  =  50  for  this  experiment. 
After  determining  the  k/s,  we  now  use  (III. 22)  to  calculate  the  C/s: 


(III. 32) 


where  Cj.  is  the  Higdon  parameter 
resulting  C/s  for  this  experiment  are: 


: 

2.2361 

H2  : 

2.0001 

2.2361 

^3  = 

2.0001 

2.0004 

2.2361 

2.0001 

2.0002 

2.0008 

2.2361 

Hb  : 

2.0001 

2.0002 

2.0004 

2.0014 

Hz  : 

2.0001 

2.0002 

2.0003 

2.0006 

H7  : 

2.0001 

2.0001 

2.0002 

2.0004 

Hs  : 

2.0001 

2.0001 

2.0002 

2.0003 

H9  : 

2.0001 

2.0001 

2.0002 

2.0002 

The  error  measurements  are  reported  in 


for  a  J^-order  boundary  condition.  The 


2.2361 

2.0021  2.2361 
2.0008  2.0030  2.2361 
2.0005  2.0011  2.0039  2.2361 
2.0004  2.0007  2.0014  2.0049  2.2361 
Figure  12  and  represent  a  significant  im¬ 


provement  over  those  reported  in  the  first  experiment.  We  now  seek  to  improve  the 
estimation  of  C/ s  using  a  method  suggested  by  Neta  and  Givoli  [Ref,  24], 


3.  Experiment  Three:  Estimating  Cj's  from  Wave 
Numbers  Obtained  using  the  Minimax  Formula 

The  third  experiment  utilizes  a  method  that  computes  C/s  from  wave  num¬ 
bers  obtained  using  the  minimax  formula  based  on  the  Chebyshev  polynomial,  an 
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Figure  12.  Experiment  2,  ||?7.h||  att=100:  Cj  from  k  Evenly  Distributed  on 
(C0=2,f=l) 


approach  adopted  by  Neta  and  Givoli  [FLef .  24].  We  consider  the  method  for  the  one¬ 
dimensional  problem  presented  by  the  current  experiment  in  which  wave  numbers  are 
given  by: 


fcj.  — 


k2 


\  2 


i  =  1, 2,  ...j  —  1, 


(III. 33) 


where  kj.  is  the  wave  number  calculated  for  a  jth-order  boundary  condition.  As 
before,  the  maximum  resolvable  wave  number  equals  50.  The  Higdon  parameters 
are  calculated  using  (III. 32).  The  J™  Higdon  parameter  Cj  is  set  equal  to  C0.  The 
resulting  Higdon  parameters  were: 


ffi  : 

2 

H2  : 

2 

2.0002 

^3  : 

2 

2.0001 

2.0007 

Hi: 

2 

2.0001 

2.0002 

: 

2 

2.0001 

2.0001 

Hq  : 

2 

2.0001 

2.0001 

H7  : 

2 

2.0001 

2.0001 

H8  : 

2 

2.0001 

2.0001 

Hg  : 

2 

2.0001 

2.0001 

2.0015 

2.0003 

2.0026 

2.0002 

2.0005 

2.0041 

2.0002 

2.0003 

2.0007 

2.0001 

2.0002 

2.0004 

2.0001 

2.0002 

2.0002 
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2.0059 

2.0009  2.0080 
2.0004  2.0012  2.0104 


The  error  nouns  reported  in  Figure  13  (left)  are  not  the  predicted  improvement  over 
the  results  reported  in  the  second  experiment.  However  with  a  minor  adjustment  to 
the  method,  we  can  bring  this  about.  Rather  than  setting  Cj  =  C0,  we  nse  (III. 29) 
and  let: 

„  _  +  P 

The  Higdon  parameters  are  now  shifted  slightly  as  follows: 


H,  : 

2.2361 

H2  : 

2.0002 

2.2361 

Hz  : 

2.0001 

2.0007 

2.2361 

tfi: 

2.0001 

2.0002 

2.0015 

2.2361 

Hz  : 

2.0001 

2.0001 

2.0003 

2.0026 

2.2361 

HQ  : 

2.0001 

2.0001 

2.0002 

2.0005 

2.0041 

2.2361 

H?  : 

2.0001 

2.0001 

2.0002 

2.0003 

2.0007 

2.0059 

2.2361 

Hz  : 

2.0001 

2.0001 

2.0001 

2.0002 

2.0004 

2.0009 

2.0080 

2.2361 

Hs  : 

2.0001 

2.0001 

2.0001 

2.0002 

2.0002 

2.0004 

2.0012 

2.0104 

The  decrease  in  error  measure  reported  in  Figure  13  (right)  indicates  that  though  the 
change  in  the  method  is  minor,  the  results  are  substantial  and  are  a  small  improvement 
over  the  results  of  the  second  experiment.  We  now  conduct  a  final  experiment  that 
will  fine  tune  the  results  found  in  this  experiment  and  minimize  the  reflected  wave. 

4.  Experiment  Four:  A  Procedure  for  Optimizing 

The  fourth  experiment  utilizes  a  procedure  for  optimihrig  Cj’s  that  can  be 
utilized  if  an  exact  solution  is  known.  It  is  computationally  intensive,  but  provides 
a  quantitative  comparison  for  the  methods  introduced  thus  far.  Imagine  the  Higdon 
parameters  Cj  to  be  controlled  by  the  dials  on  an  oscilloscope.  A  boundary  condition 
of  order  J  is  represented  by  J  dials.  The  screen  on  of  the  oscilloscope  represents 
the  finite  domain  and  displays  the  reflected  wave.  An  experimenter  adjusts  the  first 
dial  up  or  down  in  order  to  reduce  the  size  of  the  reflected  wave.  He  continues  his 
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I 

I 


Figure  13.  Left:  Experiment  3a,  ||rfc||  at  t=100  and  Cj  =  C0.  Right:  Experiment 

_  \J  o  "rntfi  /  ^  Cases,  ij.  is  Computed  from 


3b,  1 1 77 j?  1 1  at  t =100  and  Cj  = 


kfnt 


Wave  Numbers  Determined  using  the  Minimax  Formula  Based  on  the  Chebyshev 
Polynomial 


adjustments  to  the  remaining  dials  minimizing  the  reflected  wave  each  time.  After 
finishing  he  repeats  the  process  as  many  times  as  necessary  until  he  can  no  longer 
reduce  the  displayed  wave.  The  resulting  dial  settings  are  the  Cj' s  that  minimize 
the  reflection  at  the  artificial  boundary  for  a  particular  wave  action  in  the  domain  of 
interest. 

This  procedure  was  simulated  with  MATLAB  and  using  parameters  deter¬ 
mined  from  the  results  in  the  third  experiment  as  the  initial  Cj's.  The  resulting 
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Higdon  parameters  were: 


H,  : 

2.2229 

H2  : 

2.0481 

2.2359 

H2  : 

2.0170 

2.0607 

2.2361 

2.0075 

2.0260 

2.0615 

2.2361 

Hb  : 

2.0038 

2.0135 

2.0275 

2.0615 

HQ  : 

2.0031 

2.0116 

2.0206 

2.0278 

H7  : 

2.0021 

2.0075 

2.0151 

2.0276 

H8  : 

2.0012 

2.0047 

2.0094 

2.0156 

H9  : 

2.0010 

2.0037 

2.0075 

2.0106 

2.2361 

2.0616 

2.2361 

2.0608 

2.0615 

2.2361 

2.0276 

2.0615 

2.0615 

2.2361 

2.0156 

2.0276 

2.0615 

2.0616  2.2361 

The  error  measurements  reported  in  Figure  14  show  a  significant  decrease  when  com¬ 
pared  to  those  generated  by  the  modified  minimax  method.  Clearly  the  “oscilloscope” 
procedure  was  effective  in  substantially  reducing  the  magnitude  of  the  reflected  wave. 
It  is  also  interesting  to  note  that  the  maximum  Cj  for  Hz  through  H9  is  based  on 
which  is  1  in  this  experiment.  This  has  the  effect  of  allowing  the  wave  with  the 
smallest  wave  number  to  pass  through  the  boundary  unhindered.  This  is  significant 
because  this  wave  has  largest  reflection.  We  can  obtain  the  remaining  wave  num¬ 
bers  that  pass  without  reflection  at  the  boundary  using  the  above  Cj’s.  These  are 
calculated  as  follows:  _ _ 


kj.  — 


N 


C2  C2 ' 


(III. 34) 
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Using  this  formula,  the  resulting  k/s  from  the  oscilloscope  procedure  are: 


H i  : 

1.0308 

H2  : 

2.2662 

1.0004 

Hz  : 

3.8287 

2.0139 

1 

ffi: 

5.7711 

3.0908 

2.0008 

1 

Hb  : 

8.0799 

4.2890 

3.3035 

2.0008 

1 

HQ  : 

9.0027 

4.6387 

3.4745 

2.9870 

1.9992 

1 

H 7  : 

11.011 

5.7748 

4.0590 

3.0010 

2.0132 

2.0008 

1 

Hs,  : 

14.333 

7.3193 

5.1511 

4.0014 

2.9986 

2.0009 

2.0002 

1 

Hg  : 

15.918 

8.1628 

5.7489 

4.8542 

3.9942 

3.0011 

2.0004 

1.9992  1 

This  information  by  itself  is  uninteresting,  but  is  useful  when  compared  to  the  kf  s 
obtained  by  the  modified  minimax  method  which  were: 


Hi  : 

1 

H2  : 

35.355 

1 

H3  : 

46.794 

19.134 

1 

tfi: 

48.296 

35.355 

12.941 

1 

Ht  : 

49.039 

41.574 

27.779 

9.7545 

1 

Hq  : 

49.384 

44.550 

35.355 

22.700 

7.8217 

1 

H 7  : 

49.572 

46.194 

39.668 

30.438 

19.134 

6.5263 

1 

Hs  : 

49.686 

47.194 

42.336 

35.355 

26.602 

16.514 

5.5982 

1 

Hg  : 

49.759 

47.847 

44.096 

38.651 

31.720 

23.570 

14.5142 

4.9009  1 

These  results  indicate  that: 

•  It  is  more  effective  to  select  C^’s  that  correspond  to  lower  wave  numbers. 

•  The  span  of  wave  numbers  generated  by  the  minimax  method  is  much  greater 
than  the  optimized  span  indicated  by  the  “oscilloscope”  procedure. 

•  The  wave  numbers  inside  the  minimax  span  increase  too  quickly  when  com¬ 
pared  to  the  optimal  span. 
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Figure  14.  Experiment  4,  ||  at  t=100:  Cj  Optimized  using  Oscilloscope  Procedure 


In  the  next  chapter  we  explore  the  properties  of  the  Higdon  NRBC  with  respect 
to  the  two- dimension al  infinite  wave  guide. 
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IV.  DISCRETIZING  THE  SEMI-INFINITE 
WAVE  GUIDE  PROBLEM  WITH  ARTIFICIAL 

BOUNDARIES 

In  this  section  we  discretize  the  Klein-Goudon  equation.  We  also  derive  two 
forms  of  the  discrete  Higdon  NRBC  using  the  two-  and  three-point  backward- difference 
equations  for  the  first-order  derivative.  A  scheme  proposed  by  Givoli  and  Neta  [Ref. 
22]  is  used  to  simplify  this  step.  We  then  compare  the  two  Higdon  forms,  as  well  as 
previously  discussed  schemes  for  selecting  Cj,  via  numerical  example.  We  accomplish 
this  by  numerically  solving  the  semi-infinite  wave  guide  problem  (Figure  9).  A  spatial 
grid  for  j  and  y  and  a  temporal  grid  for  t  is  set  up  and  he  r/  value  at  a  point  (  jp,  yq ) 
at  time  t  =  tn  is  denoted  by: 


=*?(*P,  ?/<?,**)•  (IV. 1) 

The  rj  values  of  the  wave  guide’s  interior  points  are  determined  using  the  discrete 
Klein- Gordon  equation.  We  then  obtain  77  values  on  the  north  and  south  boundaries. 
Neumann  conditions  are  imposed  here  and  discretization  schemes  are  required.  Fi¬ 
nally,  the  rj  values  on  the  east  boundary  are  determined  using  the  discrete  form  of 
the  Higdon  NRBC. 


A.  DISCRETIZING  THE  KLEIN-GORDON  EQUATION 

The  homogeneous  linearized  shallow  water  with  zero  mean  or  Klein- Gordon 
equation  is  given  by  (11.88): 


d2V 
dt 2 


-CSV2t7  +  /2t7  =  0. 


Higdon  proved  that  in  the  context  of  the  scalar  Klein-Gordon  equation,  discrete  Hig¬ 
don  NRBC’s  are  stable  if  such  a  standard  second-order  central- difference  scheme  is 
used  to  discretize  (V.19)  [Ref.  21].  Thus  we  use  the  following  approximation  for  the 
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second  derivative: 


F"(*o)  =  lF(*o  -  h)  -  2 F(x„)  +  F(i0  +  A)]  +  0(fc!),  (IV .2) 


where  fris  the  grid  size  and  0(h2)  is  an  appioximation  euoi  of  order  k2.  We  approx- 
imate  the  Klein-Goidon  equation  by: 


1 

l W 


w, 


tj— l 
P4 


,TJ  +  1 


-  f  “  2t£?  +  ^M+l)  +  -  °- 


(IV. 3) 


Solving  fox  77p_^L  gives: 


(IV  .4) 


We  use  this  equation  to  appxoxixnate  the  intexiox  points  of  the  doxnain.  The  appxood- 
mation  exxox  is  of  oxdex  0(  Ax2,  Ay2,  At2).  Note  that  the  schexne  is  two  level  in  time 
and  thus  one  xe quires  an  additional  starting  value.  This  can  he  accomplished  using 
the  Matsuno  scheme  [Ref.  33], 


B.  DISCRETIZING  THE  NORTH  AND  SOUTH 
BOUNDARIES 


If  the  conditions  on  the  north  and  south  boundaries  contain  derivatives,  then 
they  too  must  be  discretized.  Fox  example,  consider  the  Neumann  type  condition: 

The  south  boundary  is  discretized  by  using  a  three-point  forward- difference  formula 
for  F'(x0): 

F'(x0)  =  -!-[-3F(x0)  +  4F(x0  +  h)  -  F(x0  +  2 h)  +  0(h2).  (IV. 5) 


drj 
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Thus  the  south  "boundary  condition  is  approximated  by: 


%=^y  |-3<»  +  _  +  °( V)  =  °>  (™.S) 

where  5  is  the  y- index  marking  the  south  boundary,  A y  is  the  size  of  the  y-grid,  and 
0(  Ay2)  is  an  approximation  error  of  order  Ay2.  We  now  approximate  rj  on  the  south 
boundary: 

^.s  =  |<s+i  -  ^.s+2  •  (IV  ■7) 

Similarly  on  the  north  boundary,  three-point  backward- difference  formula  for  F'(x0), 
an  approximation  of  order  Ay2  for  77  is: 

Vp.n  =  ^p-Jv-i  —  2^p-Jv-2  >  (IV .8) 

where  N  is  the  y- index  marking  the  north  boundary. 


C.  THE  DISCRETE  FORM  OF  THE  HIGDON  NRBC 

The  J*h- order  Higdon  NRBC  Hj  (III .4)  is  the  product  of  J  operators: 


n(|+<^=°- 


;=1 


dx 


(IV. 9) 


To  discretize  this  equation,  use  the  two- point  backward- difference  formula  for  the  first 
derivative: 

F'(xo)  =  ±  [F(x 0)  -  F(x 0  -  h)]  +  0(h),  (IV. 10) 

where  0(h)  is  an  approximation  error  of  order  h.  The  discrete  form  of  (IV.  9  )  becomes: 

J  r  /y(x,y,t)-y(x,y,t-  At)\ 


n 

?=!•  L 


+  C) 


At 


'y(x>  y,t)-y(x  -  Ai,y,t)’ 


(iv.li) 


=  0. 


The  approximation  error  for  this  equation  is  of  order  0(Ax,  At). 
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Two  special  operators  are  defined  to  further  facilitate  the  discussion.  The 
temporal  shift  operator,  denoted  "by  Sf,  affects  a  grid  point  as  follows: 


Vlq  =  7Jp4L  =  y(*P’  «j.  -  At)- 


(IV. 12) 


Similarly,  the  spatial  shift  operator  in  the  x- direction,  denoted  by  S  ,  affects  a  grid 


point  as  follows: 


s*y?4  =  y£-u  =  y(xP  - 


(IV. 13) 


The  identity  operator,  denoted  by  I,  has  no  effect  when  operating  on  the  grid  point, 
e.g.  Using  these  operators,  the  discrete  form  of  the  Higdon  NRBC  is: 


TT  [7~S'~  i  c  (T~S* 

At  Ax 


n h  = 


(IV. 14) 


where  the  index  E  corresponds  to  the  location  of  the  Higdon  NRBC  at  r  =  xE,  By 
collecting  the  operator  terms  I,  Sf,  and  Sf,  (IV.  14)  can  be  rewritten  as: 

n  [(i  +  <^) 1  -  sr  -  (c,^)  s;] nl,  =  0-  (iv .15) 


Making  the  following  substitutions: 


At  At 

(j-i  =  1  (Sj ~r  ,  d,  =  — 1,  and  €«  =  — , 

J  Ax  J  3  JAx 


(IV. 16) 


allows  us  to  rewrite  (IV. 15)  as: 


j 

0  I  4-  djSt  ejS^  'j  yE_q  =  0. 


(IV. 17) 


Now  consider  a  second  order  Higdon  NRBC  H2,  which  by  using  (IV.  17  )  and 


expanding  can  be  written  as: 

-\-d1d2ISf  — )- O-l €2lSx  -\-d]_d2lSi  -\-did2St  2 

-| -die2SfSf  -\-eia2ISf  -fei  d2Sf  Sf  +ei  e2Sf2]7jEjq  =  0. 

In  this  form,  Hj  is  represented  as  the  summation  of  3J  terms: 

3  J  — L 

Y. ]  4^4 ?  = 

m=0 


(IV  .18) 


(IV  .19) 
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where  Am  is  a  product  of  a j,  dj  and/ox  ej,  and  Pm  is  made  up  of  a  combination  of 
operators  I,  Sf  and/or  S~ .  The  first  term  of  (IV. 17)  is: 

(IV.20) 


By  letting: 

J  3J-L 


a  =  n 

and  Z*  =  Y.  A^P^Ves}  » 

(IV. 21) 

2=1 

m= L 

we  can  rewrite  (IV. 19)  as: 

A  +  Z*Ve4  ~  0- 

(IV. 22) 

Since  Ac  4-  0,  this  implies: 

•h 

II 

1 

Sb.  *2 

°  ft 

•tb 

(IV. 23) 

The  problem  now  is  to  evaluate  Z*. 

We  now  employ  the  scheme  devised  by  Givoli  and  Neta  [Ref,  22],  Consider 
the  coefficient  Am  and  the  operator  Pm .  We  rewrite  the  index  m  in  base  3  and  refer 
to  it  as  rn(3).  Table  II  summarizes  An(j,  a-nd  and  their  corresponding  values  for 
the  discrete  form  of  Hi  described  by  (IV. 18).  Inspection  of  the  table  reveals  a  useful 


Table  II.  H2  Values  for  Am^  and  Pm(3;] 


A  - 1  A»ta)  =  ala2 

Pq  — 1 •  Poo(3,  —  J2 

'* 

II 

n' 

<s 

Ai  —*■  Aif3,  =  aidi 

Pl  — >  Pol(3,  =  ISt 

P°^Ve,  =  Ve7 

Al  —*  ^4o2(3)  = 

P2  — *•  Po2(3,  =  IS~ 

Pni3)VE.Q  =  Ve-Ij} 

Pz  -»•  P10,3,  =  St~[ 

S3 

tj 

ft 

A 

11 

ft 

•k,  1 

1- 

A  Alfa,  =  dldi 

Pu-tnVE.q  -  Ve/ 

A$  —*■  A2(a,  =  die2 

P&  — 5 ‘  -Pl2(3)  — 

P^mVE-Q  —  Ve-Ijj 

A  ~ ^  Ao(a,  =  ela2 

P&  P 20(3,  =  S?  I 

1 

II 

Or 

I? 

<s 

A  — Alfa,  =  e^d'i 

P7  -*■  ^21,3,  =  SxSf 

— 1  •— 1 

II 

Or 

5 

cC 

Ag  —*■  ^22,3,  =  ^1^2 

m.  bi.mjhm 

P^tn^E-a  -  Ve-24 

pattern  that  allows  us  to  quickly  generate  the  discrete  form  of  Hj  for  any  order  J. 
With  regards  to  the  index  m(3) ,  a  digit  value  of  0,  1,  ox  2  implies  that  the  coefficients 
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Table  III.  Constructing  Pm  and  A m  from  the  Base  Thiee  Index  0221(3) 


Digit  Position 

1 

2 

3 

4 

Digit  Value 

2 

2 

1 

Use  Digit  Value  to 

Select  Constant 

(Zj 

Use  Digit  Position  to 

Set  Constant’s  Subscript 

dy 

<?2 

Use  Digit  Value  to 

Select  Operator 

I 

C-l 

St~l 

dj,  dj,  or  ej  respectively  are  part  of  the  product  that  make  up  An(3).  The  position  of 
the  digit  determines  the  index  j  for  each  term,  where  we  consider  the  left  most  digit 
to  be  in  position  1.  Similarly,  the  digits  0,  1,  and  2  imply  the  application  of  the  I, 
Sf,  or  S~  operator  respectively  in  the  construction  of  •  For  example,  H ^  has 
81  terms.  The  index  for  the  25th  term  has  a  base  3  equivalent  equal  to  0 22 1  j  (note 
that  J  digits  are  always  used  for  Hj).  Table  III  summarizes  the  application  of  the 
base  three  index  0221(3j  to  construct  Am  and  Pm.  Therefore: 

^25  ~ 1  A)22L(a)  =  ale2e2^-l  and  P26  — ►  P 5221(3j  =  > 

and: 

We  can  further  generalize  the  operator  equation  as  follows: 

=  iSX,  ,  (IV  .24) 

where  d  and  b  are  the  number  of  times  that  the  digits  1  and  2  appear  respectively  in 
™(3)- 

Thus  far  we  have  indicated  that  this  scheme  will  generate  3J  terms,  but  this 
number  is  reduced  considerably  when  one  combines  terms  containing  the  same  opera¬ 
tor.  For  example,  in  Table  II  the  operators  ISf ,  IS~,  and  5“  Sf  appear  twice.  When 
combined,  the  number  of  terms  for  the  discretized  form  of  P2  Is  reduced  by  three,  In 
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general,  terms  will  contain  operators  of  the  form  Sx  aSt  b  where  0  <  a  4-  b  <  J.  The 
number  of  integer  combinations  of  a  and  b  that  satisfy  the  inequality  is: 

( J  +  2)(  J  +  1) 


j+i 
3=1 


2 


(IV. 25) 


and  represents  the  number  of  different  operator  possibilities  for  Hj,  This  is  significant, 
because  it  greatly  reduces  storage  requirements  as  well  as  processing  time  when  using 
higher  order  Higdon  NRBC’s.  For  example,  is  fully  discretized  using  55  vice  19,683 
(39)  terms. 


D.  IMPROVED  FORM  FOR  THE  DISCRETE  HIGDON 
NRBC 

In  the  previous  section,  we  formulated  a  discrete  form  of  the  Higdon  NRBC 
with  an  approximation  error  of  order  0(h).  The  approximation  error  for  the  dis¬ 
cretized  Klein- Gordon  equation  and  for  the  north/south  boundaries  were  of  order 
0(h2).  This  suggests  that  we  might  improve  our  scheme  for  discretizing  the  Higdon 
NRBC  by  using  the  three-point  backward- difference  formula  which  is  given  by: 


F'(x0)  =  —  [3F(  j0)  -  4F(j0  -  h)  +  F(j0  -  2h)\  +  0(h2), 

Ih 


(IV. 26) 

Using  this  in  operator  notation,  the  improved  discrete  form  of  the  Higdon  NRBC  is: 

>£,  =  0-  (IV  .27) 


J  r3/  -  4Sf  +  Sf^  +  /3£  -  45“  +  S-!  ’ 


n 

3=1 


At  '  3  V  Ax 

Collecting  the  operator  terms  J,  Sf,  and  S~,  (IV. 27)  can  be  rewritten  as: 

=  t™'28) 


JL  [/  A t\  4  1  ,  4  At  1  At 

n  K1  +CiA;)  1  ~  Ist  +  3S,_  “  3CjAjS'  +  %C>  a7S>  . 

Making  the  following  substitutions: 

.  _  At  4  1  4  At  1  At 

*3-1  +  ^^,  bj  3 ,  Cj-3,  dj  3 Oj  a ^ j  and  ^-3^^, 

allows  us  to  rewrite  (IV. 2 8)  as: 

/ 

n  (a3^  +  hsr  +  c35r2  +  2)  ifEjq  =  o. 

3=1 


(IV. 29) 


(IV. 30) 
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In  the  expanded  form,  Hj  is  represented  as  the  summation  of  5J  terms: 


5  J  — L 


(IV. 31) 


m=0 


where  is  a  product  of  dj,  bj,  Cj,  dj,  and/or  ej.  Similarly  Pm  is  made  up  of  a 
combination  of  operators  I,  Sf,  Sf2,  S~  and/or  S~2.  By  letting: 


J  bJ-l 

Aq  =  J  J  dj  and  Z  =  )  '  AmP j 

^=1  m=l 

and  noting  that  jlo  ^  0  we  rewrite  (IV. 31)  as: 


Ve*  = 


Aq 


(IV.  32) 


(IV. 33) 


Using  a  procedure  similar  to  that  employed  by  Givoli  and  Neta  [Ref.  22],  we  can 
easily  sift  through  the  algebraic  complexities  of  (IV.  33). 

Consider  again  the  coefficient  Am  and  the  operator  Pm .  We  rewrite  the  index 
m  in  base  5  which  we  will  refer  to  as  .  With  regards  to  this  index,  a  digit  value  of 
0,  1,  2,  3,  and  4  implies  that  the  coefficients  dj,  bj,  Cj,  dj,  or  ej  respectively  are  part 
of  the  product  that  malces  up  ■  The  position  of  the  digit  determines  the  index  j 
for  each  term,  where  the  left  most  digit  is  in  position  1.  Similarly,  the  digits  0,  1,  2, 
3,  and  4  imply  the  application  of  the  I,  Sf,  Sf2,  S~,  or  S~2  operator  respectively  in 
the  construction  of  Pm(4) ,  For  example,  Ht,  has  3125  terms.  The  index  of  the  1454th 
term  has  a  base  5  equivalent  of  21304(5)  •  Therefore: 

P LJ54.  =  , 

and: 

^4uSJ.Pu64.T/£:4  = 

We  have  indicated  that  this  scheme  will  generate  5J  terms,  but  this  number 
is  reduced  considerably  by  combining  terms  with  the  same  operator.  Terms  will 
contain  operators  of  the  form  IS~aSfb  where  0  <  a-|-f>  <  2J.  The  number  of  integer 
combinations  of  a  and  b  must  be  less  than; 

(2J+2)(2J+1) 


2J+1 
j=l 


2 


(IV. 34) 
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We  note  that  certain  opeiatoi  combinations  are  not  possible.  For  example,  consider 
Hz.  The  operator  SfbS~  is  impossible  because  all  three  operators  are  required  to 
construct  the  Sfb  term  leaving  no  possibility  for  it  to  coexist  with  the  S~  term. 
Taking  into  consideration  these  operators,  the  number  of  terms  needed  to  describe 


Hy. 


(2J+1)2  +  1 

2 


Hence  Hg  is  fully  discretized  using  181  vice  1,953,125  (59)  terms. 


(IV.35) 


E.  NUMERICAL  EXAMPLE:  EMPLOYING  HIGDON  NRBC’S 
FOR  A  SINGLE  WAVE  PULSE  AT  EV  IN  A  SEMI¬ 
INFINITE  WAVE  GUIDE 


We  now  return  to  the  semi- infinite  wave  guide  depicted  in  Figure  9  of  Sec¬ 
tion  A.  We  let  the  wave  guide  width  6=5  and  depth  ©  =  .1.  The  medium  is 
dispersive  with  /  =  .5  and  a  gravitation  acceleration  of  g  =  10  is  used.  The  initial 
values  are  zero  everywhere.  The  boundary  function  on  the  west  boundary  fV  is: 


*7»-(y>*)  = 


|  .005©  cos 

n  .  ' 

2 r{V  ~  Vo\ 

1 0 

(IV. 36) 


0  <  t  <  fo 

otherwise  , 

where  r  and  t0  are  the  wave  pulse  center,  radius  and  time  duration  respectively. 

The  parameter  values  are  set  at  tfo  =  2.5,  r  =  1.5,  and  f0  =  0.5. 

An  artificial  boundary  G  is  imposed  at  x  =  5,  defining  the  computational 

domain  ft  as  a  5  X  5  square.  A  mesh  of  20  X  20  is  used  in  ft.  The  two-point  backward 

difference  method  as  described  in  Section  IV.  Equation  IV.  23  is  used  to  estimate  the 

Higdon  boundary  values.  The  extended  domain  V  for  the  reference  solution  77  is 

a  15  x  5  rectangle  with  a  60  x  20  mesh,  An  artificial  boundary  is  imposed  on  V  at 

x  =  15.  Any  spurious  reflections  resulting  from  this  boundary  will  not  affect  17 

25 

if  the  run  time  is  less  than  — — .  This  the  time  it  takes  for  a  wave  to  travel  from 

Q) 

x  =  0  (r^)  to  x  =  15  (rv  for  V)  plus  the  time  it  takes  for  any  reflection  to  travel 

&n 

back  to  x  =  5  (TH,  fo*  ^).  On  the  north  and  south  boundaries  we  impose  —  =  0. 
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This  “hard  walT  condition,  which  causes  waves  to  "bounce  off  the  north  and  south 
boundaries,  will  result  in  additional  geometric  dispersion. 

Two  cases  with  artificial  boundaries,  17^,  and  17^,,  are  computed  and  jux¬ 
taposed  to  7 7  .  For  77mm1,  an  NRBC  with  J  =  4  is  constructed  with  parameters 

Cj  =  {Co,  Co,  Co,  Co} .  Here  Co  =  y/ g@  =  1  is  the  minimum  wave  speed.  An  error 
measurement  error  ||g’(f)||r,  on  the  truncated  domain  ft  is  obtained  by  using: 


JV, 

=  hy. 

j=l  j=  1 


NXN« 


(rv.37) 


where  Nx  and  Ny  are  the  number  of  grid-points  in  the  x-  and  y-  directions  as  de¬ 
termined  by  the  grid  spacing.  For  77^,,  an  NRBC  with  J  =  1  and  Cj  =  {C0}  is 
constructed  and  its  numerical  solution  is  compared  to  77^  to  obtain  a  second  er¬ 
ror  measurement.  We  rate  the  effectiveness  of  NRBC’s  by  comparing  the  values  of 
||e'(t)||r, .  Note  that  ||e(t)||a  is  a  function  of  time  and  will  start  to  increase  as  the  wave 
impinges  on  the  artificial  boundary. 

Figures  15  through  19  show  the  solutions  for  77^,  77^^  and  77^,  at  times 
i=l,  4,  5,  8  and  10  seconds.  The  top-left  and  -right  plots  depict  77^  on  the  truncated 
domain  ft  and  extended  domain  V  respectively  (note  that,  though  the  extended 
domain  is  continuous,  it  is  separated  in  the  figure  so  that  77  may  be  better  contrasted 
with  T7ca-1  and  77...., ),  The  middle-  and  bottom-left  plots  correspond  to  77mm1  and 
respectively,  Two  graphs  on  the  center-  and  bottom-right  present  ||e(t)||n  for  77mm1 
and  77mmJ  .  Note  that  “HNRBC-2DR-1S-1L-U0-V0-T0T  appears  in  the  caption.  This 
shorthand  will  be  used  throughout  the  paper  to  identify  to  identify  the  figures.  It  is 
defined  as  follows: 


•  HNRBC:  Higdon  non-reflecting  boundary  condition, 

•  tiDR:  ti- dimensional  rectangular  domain, 

•  tiS:  HNRBC  applied  to  n  sides  of  domain, 

•  nL:  Ti-1  aver  problem, 
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•  Umpn:  x-directdon  advection  U  =-.n  (note:  m  denotes  a  minus  sign  and  p  a 
decimal  point) 

•  Vpn:  ^-direction  advectaon  V  =.n. 

•  T n:  solution  at  t  —n. 

Therefore  “HNRBC-2DR-1S-1L-U0-V0-T01”  lefeis  to  the  solution  at  t  =  1  of  a  prob¬ 
lem  in  which  a  Higdon  NRBC  as  applied  to  one  side  of  a  two-dimensional  rectangular 
domain  consisting  on  a  single  layer  with  no  advectaon,  We  will  introduce  multi-layer 
stratification  and  advection  models  in  later  chapters. 


ti  ,an  Enisiar  Omar  ,  _a 

x  10 


Figure  15.  HNRBC-2DR-1S-1L-U0-VO-T01  with  Simple  Pulse  Boundary  Condition 


At  t  =  1  (Figure  15)  the  wave  packet  is  still  close  to  r«.-  and  no  spurious 
reflections  have  occured.  The  plots  for  77  , 77^^  and  77^^,  are  identical  and  ||e(t)||n 

for  77  ,  and  77  ,  is  essentially  0. 

luisl  lcai«2 

At  t  =  4  (Figure  16)  the  leading  edge  of  the  wave  packet  reaches  B.  A  slight 
spurious  reflection  is  measured  for  T7caiil  and  77^^,  but  overall  the  three  solutions  are 
still  very  similar.  At  t  =  5  (Figure  17)  the  wave  packet  crosses  B.  The  77 _ ,  and 
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Figure  16,  HNRBC-2DR-1S-1L.-U0-V0-T04  with  Simple  Pulse  Boundary  Condition 


V^.^1  P^ts  are  still  indistinguishable  from  tj  ,  At  t  =  8  and  10  (Figures  18  and  19), 
most  of  the  wave  packet  has  left  and  is  now  visible  in  the  extended  domain  I?.  The 
T7m-1  solution  exhibits  wave  traces  similar  to  those  in  77  ,  whereas  the  17^,  solution 

does  not,  The  difference  in  scale  for  ||  )  || n  reveals  an  improvement  of  one  order 
of  magnitude  for  This  is  quantitatively  significant,  but  the  qualitative  results 

are  of  greater  note,  For  t^m1,  the  wave’s  energy  has  passed  through  B  relatively 
unimpeded.  On  the  other  hand,  reveals  visible  spurious  reflections  resulting  in 
a  backwards- moving  wave  that  pollutes  the  computational  domain 

This  example  illustrates  the  robustness  of  the  Higdon  NRBC,  The  perfor¬ 
mance  of  the  NRBC  is  improved  by  simply  increasing  the  order  J.  We  now  use  this 
example  to  measure  any  improvement  in  the  solution  using  the  two-point  Higdon 
NRBC  approximation  derived  in  Section  IV. C  versus  the  three-point  Higdon  NRBC 
approximation  derived  in  Section  IV, D,  In  both  cases  Cj  —  (Co,  Co,  C0,  Co}  where 
Co  =  1,  The  results  in  Figure  20  report  a  solution  improvement  of  an  order  of  about 
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Figure  17.  HNRBC-2DR-1S-1L-U0-V0-T05  with  Simple  Pulse  Boundary  Condition 


1/3  when  three- point  approximations  are  used.  We  also  use  the  current  example’s 
parameters  to  compare  results  when  the  Cj  are  chosen  as  a  preprocess.  An  algo¬ 
rithm  using  the  symmetric  minimax  formula  (based  on  the  Chebyshev  polynomial) 
proposed  by  Sommeijer  et  al,  [Ref.  32],  In  Case  1,  the  Higdon  NRBC  has  an  order 
J  =  4  with  Cj  =  {Co,  Co,  Co,  Co}  while  in  Case  2,  J  =  4  and  the  Cj’s  are  chosen 
as  a  preprocess.  In  both  cases  the  2nd  order  approximation  for  the  Higdon  NRBC  is 
used.  The  expected  improvement  is  not  revealed  by  Figure  21.  This  discrepancy  will 
be  checked  in  a  follow-up  example. 
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Figure  18.  HNRBC-2DR-1S-1L-U0-V0-T08  with  Simple  Pulse  Boundary  Condition 
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Figure  19.  HNRBC-2DR-1S-1L-U0-V0-T10  with  Simple  Pulse  Boundary  Condition 
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Figure  20.  ||e’(t)||n  Plots  foi  2-Point  vs.  3-Point  Higdon  NRBC  Approximations 
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Figure  21.  ||e(i)||n  Plots  for  Cj  =  {C0,  C0,  C0,  C0}  vs.  Cj  Selected  using  a  the 

Symmetric  Minimax  Formula  with  3- Point  Higdon  NRBC  Approximations 
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F.  NUMERICAL  EXAMPLE:  SEMI-INFINITE  WAVE  GUID  E 
WITH  A  CONTINUOUS  WAVE  INTRODUCED  AT  Fw 

The  set  up  foi  this  example  is  the  same  as  the  previous  example  with  the 
following  exceptions:  Instead  of  a  single  pulse,  the  boundary  function  on  IV  is  con¬ 
tinuous,  linear  combination  of  three  waves  with  the  form: 

V (y,  *)  cos  (“ipO  CC6(^x  -  ^rnt),  (IV. 38) 

m=L  \  0  / 

whose  parameters  are  selected  a  priori  to  be: 


Am  =  .01,  .01,  .01: 

n™  =  1,  2,  2: 

ti>„ ,  =  .81,  1.37,  1.68. 

To  satisfy  the  Klein- Gordon  equation,  the  dispersion  relation: 

<  =  Cl  (ki  4^4  f. 

must  be  invoked  to  determine  fc^: 

,  _  K  -  P  <*2 

^  V  Cl  62  ' 

Thus  using  the  parameters  given  for  (IV.  38),  is: 

=  .11,  .22,  1.00. 


(IV. 39) 


(IV  .40) 


Two  cases,  and  are  again  juxtaposed  to  77  In  Figure  22.  For  77. _ , . 

an  NRBC  with  J  =  5  and  parameters  Cj  =  {C0,  C0,  C0,  C0,  C0}  is  used.  For  77...., , 
J  =  2  and  Cj  =  (C0,  C0}  is  used.  An  error  measurement  || e’(t)  || n  given  by  (IV. 37)  is 
determined.  As  expected,  the  higher- order  Higdon  NRBC  results  in  asmaller  ||  ^(t )  ||  n 
and  therefore  generates  less  spurious  reflection  at  B. 

We  now  examine  previous  suppositions  considered  in  Section  III.C  concerning 
the  selection  of  the  Higdon  parameters  Cj. 
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Figure  22.  HNRBC-2DR-1S-1L-U0-V0-T10  with  Continuous  Wave  Boundary  Condi¬ 
tion 


•  In  Figure  23  we  compare  Cj  =  {1, 1, 1, 1, 1}  vs.  a  5th- order  Higdon  NRBC 
where  the  C/s  are  automatically  selected  using  the  minimax  formula  based 
on  Chebvshev  polynomials.  The  result  favors  the  former,  discounting  previ¬ 
ous  supposition  that  auto-selection  of  C^s  using  this  scheme  might  yield  an 
improvement. 

•  In  Figure  24  we  compare  Cj  =  {1, 1, 1, 1, 1}  to  a  5th-order  NRBC  where  the 
Cj’s  are  evenly  distributed  from  C0  to  yjcl  4-  /2.  Any  improvement  indicated 
by  using  the  latter  scheme  is  not  significant. 

•  In  Figure  25  we  compare  Cj  =  {1, 1, 1, 1, 1}  to  a  set  of  five  Cj's  that  are  evenly 

.  7T 

distributed  from  Co  to  the  maximum  resolvable  wave  number  ; — . 

5Aj 

Again  the  simpler  scheme  for  choosing  Cj  s  is  numerically  superior, 

Based  on  results  in  this  and  the  previous  section  we  conclude  that  spurious  reflections 
at  the  Higdon  boundary  are  reduced  by  using: 


•  Cj  s  equal  to  {C0,  ...C0}. 

•  Higher  order  Higdon  NRBC’s. 
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•  three-point  approximations  for  the  derivatives  in  the  Higdon  NRBC. 


*10^ 


Figure  23.  ||e(f)||n  for  Cj  =  {1, 1, 1, 1, 1}  vs.  Cj  Auto-Selected  Using  Minimax  For¬ 
mula  Based  on  Chebyshev  Polynomials. 

We  make  one  final  comparison  before  closing.  In  Figure  26  a  Higdon  NRBC 
with  J  =  5  that  uses  a  two-point  approximation  is  juxtaposed  to  a  Higdon  NRBC  with 
J  =  4  that  uses  a  three-point  approximation.  We  see  that  the  results  are  quite  close. 
However  differences  in  computational  effort  is  significant.  From  Section  IV.  C  we  know 
that  to  discretize  a  two-point  Higdon  NRBC  with  J  =  5  required  the  generation  of 
125  (36)  terms  that  were  subsequently  reduced  to  21  terms  when  redundancies  are 
combined  (IV. 25).  From  Section  IV. D  we  determined  that  to  discretize  a  three- 
point  Higdon  NRBC  with  J  =  4  required  the  generation  of  625  (511)  terms  that 
were  subsequently  reduced  to  41  terms  (IV. 35),  Since  the  discrete  Higdon  NRBC 
formula  must  be  applied  every  time  an  artificial  boundary  grid  paint  is  encountered, 
it  is  apparent  from  this  comparison  that  it  is  more  efficient  to  use  the  two-point 
approximation  with  a  higher  Higdon  order  J. 
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Figure  24.  ||e(t)||n  for  Cj  =  {1, 1, 1, 1, 1}  vs.  Cj  Distributed  Evenly  from  C0  to 


Figure  25.  ||e(f)||n  foi  Cj  =  {1, 1, 1, 1, 1}  vs.  Cj  Distributed  Evenly  from  C0  to 

k™**  =  bKx 
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Figure  26.  ||e(t)||n  for  Cj  =  {1, 1, 1, 1, 1}  Using  a  Two-Point  HNRBC  Approximation 
vs.  Cj  =  {1, 1, 1, 1}  Using  a  Three-Point  Higdon  Approximation 
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V.  A  N-LAYER  STRATIFIED  DISPERSIVE 

WAVE  MODEL 


We  now  lift  the  homogeneous  fluid  assumption  and  develop  equations  to  model 
geophysical  flow  in  a  stratified  medium.  A  suitable  medium  for  a  geophysical  dynam¬ 
ics  is  the  open  ocean  where  fluid  density  p  is  affected  by  salinity  and  temperature. 
Salinity  changes  are  slight  in  this  environment,  and  temperature  remains  relatively 
constant  in  the  horizontal  directions.  However,  temperature  does  change  significantly 
in  the  vertical  direction.  Therefore  we  approximate  p  to  be  a  function  of  z  only.  Since 
the  fluid  is  assumed  to  be  incompressible,  p  does  not  vary  with  pressure  p. 

The  linearized  shallow  water  equations  (11.78  and  11.80)  were  derived  in  part 
from  the  continuity  equation  for  homogenous,  incompressible  fluids  (II. 6): 


du  dv  dw 

^  +  ai  +  si =  °- 


It  was  shown  by  (11.59)  that  u  and  r  are  independent  of  z  for  a  constant  density 
fluid  enabling  us  to  uncouple  w  from  u  and  r  in  (II  .6),  which  after  integrating  yielded 
(11.61): 

w(x,y,z,t)  =  -z 

This  critical  step  in  the  derivation  is  no  longer  possible  when  we  assume  that  p  is 
dependent  on  z. 

We  extricate  ourselves  from  this  conundrum  by  developing  a  layered  shallow 
water  approximation  where  p  is  constant  in  each  layer  (Fig.  27).  Here  it  is  assumed 
that  the  fluid  is  still  incompressible  and  that  density  pi  is  constant  in  each  layer  Li, 
but  varies  in  the  different  layers.  In  order  for  this  stratification  scheme  to  be  stable, 
pi  must  be  monotonically  increasing  downward  [Rjef.  28].  Additionally  we  assume 
that  there  is  no  fluid  mixing  between  layers. 


y,t)  + 


dx 


$y 


-1-  w(x,y,  t). 
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Figure  27,  AT- Layer  Shall ow  Water  Model 


Referring  to  the  AT- layer  shallow  water  model,  the  pressure  p{  at  any  point  in 
Lj  is  determined  from  hydrostatic  principles: 


't-i 


JV 


Pt  =  Po  +  g\^,Pjhi-\-fn'£'(hj-z)\f  (V.l) 

\;=1  j=i  } 

where  Fo  is  a  constant  ambient  pressure  at  the  surface  ho  and  N  is  the  total  number 
of  layers  in  the  model.  In  (V.l),  the  first  summation  term  is  the  contribution  to  pi 
from  the  layers  above  Li.  The  second  summation  term  is  the  contribution  to  from 
the  liquid  column  in  L, .  We  use  (11.80)  and  (V.l)  to  obtain  the  horizontal  momentum 
equations  in  L*: 


68 


Derivation  of  the  vertical  momentum  equation  in  L*  is  more  complex.  Since 


pi  is  constant  in  Li,  we  can  uncouple  the  continuity  equation  for  L*: 


Wi(x,y,  z,i)  = 


dui(jc,y,t)  dvi(x,y,t)' 


-\-Wi(x,y,t), 


(V.3) 


dx  dy 

where  Wi  is  a  vertical  velocity  component  in  L*.  For  brevity  we  drop  dependent 
variables  from  subsequent  expressions.  At  the  interface  between  Lj_i  and  Li  where 

JY 

z  =  hj,  the  vertical  speed  component  Wi  (see  similar  discussion  in  Section  II. B. 6) 

j=i 

^=4x>)+*4i>.+'4x>i-  (v.4) 


atfs 

j=t 


This  implies  that: 


w4 


9  " 


tdxft 


N 


dy  j=i 


N 


^  H  hi  +  57  E  hi )  +  ar  1 v*  ]  • 


j=t 


dx 


;=t 


ay 


(V.5) 


Similarly,  at  the  interface  between  L{  and  Li+l  where  z  = 
component  Wi  is: 


iY 

V  hj,  the  vertical  speed 

J=t+L 


d  N  d  N  d  N 

w< = ^  (v'6) 

which  implies  that: 

a/JV\a/JV\ 

(V7) 

Since  ihj  is  independent  of  z,  (V.5)  and  (V.7)  must  be  equal.  Therefore: 

(V.8) 

Equation  (V.8)  is  the  vertical  momentum  equation  for  Li.  Together  with  (V.2)  this 
completes  the  description  of  the  fluid  motion  inside  of  the  ith- layer  L{. 

Before  considering  a  numerical  solution  we  linearize  the  governing  equations 
for  each  layer.  We  assume  that  the  7ih  r,,  and  h^  are  dominated  by  constant  terms 
U i,  V]  and  ©j.  Superimposed  on  these  are  small  variations  v*,  v*,  and  of  0(8),  i.e . : 

Ui  =  Ui-\-  ti*,  1'i  =  Vi  4-  r*,  and  hi  =  0*  4-  rji.  (V.9) 


d  d  d 

at11  + 
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Substituting  these  in  (V.2)  and  (V.8)  and  neglecting  terms  of  0(£2)  yields: 


dt 


( t— L 


f- + ■ u<£' + +-» = - (§ t ) .  (v'10) 


^ir/iivi  +  e  /&£  at£\ 

at  +  ^  ar  +  Vt  dy  +  &i{  dz  +  dy  )  ~  °’ 

If  we  assume  that  there  is  no  advection  (e.g.  Ui  =  V]  =  0),  then  (V.10)  reduces  to: 


(V.11) 


Equation  (V.ll)  is  the  linearized  shallow  water  equation  for  the  stratified  TV-layer 
model  with  no  advection, 

A.  REDUCING  THE  STRATIFIED  TV-LAYER  MODEL 
TO  THE  KLEIN-GORDON  FORM 

In  order  to  combine  linearized  components  of  the  stratified  TV-laver  model 
(V.ll)  to  a  Klein- Gordon  form,  we  use  a  step-by-step  derivation  that  is  similar  to 
that  used  in  Section  II, B, 8. 

Step  1:  Perform  the  following  substitutions  in  (V.ll): 

if*  =  and  v{  =  r*©*. 
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Recalling  that  0*  is  independent  of  time,  the  Jesuit  is: 


j=L  pi  Bx  dx 


)=t 


±1  a  dy  j:,dy 


(V.12) 


(V.13) 


3m  dUi  dii 

—  4-  — 1  -I-  — -  =  0. 
Bt  dx  dy 


(V.14) 

Step  2:  Let  /  be  constant.  Take  the  partial  derivative  of  (V.12)  with  respect  to  x: 


d  f  Bdi 


dv{ 


'i-l 


5 


ft  ^ 


dx 


ft  5*2  +  £ 

and  the  partial  derivative  of  (V.13)  with  respect  to  y : 


dx2 


.  du, 


i-l 


Pjd%  , 


J=Ift 


;=! 


a?/2 


Add  the  resulting  equations: 


s(^8M$bt)-^(gSy£ 

Step  3:  Take  the  partial  derivative  of  (V.12)  with  respect  to  y\ 


ft 


(V.15) 


d  ( Biii 


,Bvi 


'i-l 


"E77  [  Cu  _  f  Cl..  —  _ ■ 


ft  i  v 


By  \Bt 


dy 


By  l  j=1  pi  Bx 


+  E^ 


)=t 


and  the  partial  derivative  of  (V.13)  with  respect  to  x: 


B  ( Bii 


.  Biii 


(i-l 


=  1^1 +/£—**=  IE 


Bx  V  Bt 


Pi  drjj  ^  ^  Brij 


dx  L=1  pi  dy  j=i  By 


Subtract  the  resulting  equations: 


B  ^  Bui  Bii  ^  ^  l'  Bii  ^  Biii  ^ 


Bt  \  By  Bx 


\By  ^  Bx 


0. 


(V.16) 


Step  4:  Take  the  partial  derivative  of  (V.15)  with  respect  to  t  and  rearrange  the 
resulting  terms: 


Bt  \Bx  By 


U= l  ft 


J=! 
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Multiply  (V.16)  by  /  and  rearrange  the  resulting  terms: 


d  f  dr  j 
fat  [dx 


dvi  dUi  \ 

¥+&)' 


Add  the  resulting  equations: 


\dt*  '  J  )\dx  '  dy)  '  *  %at 
Step  5:  Using  (V.14)  rewrite  (V.17)  as: 


j=t 


d_ 

dt 


f*+/»w-*,W ga,+E* 


J=! 


=  0. 


Step  0:  Integrating  (V.18)  with  respect  to  t  yields: 


(V.17) 


(VIS) 


(V.19) 


where  Si{x,y )  is  an  arbitrary  function  that  results  from  integration.  Equation  V.19, 
resembles  the  Klein- Gordon  form  reported  in  (11.88),  and  is  a  restatement  of  the  set 
of  linearized  (about  a  zero  mean)  shallow  water  equations  (V.ll)  generated  for  the 
N-layer  stratification  model. 


B.  DISCRETIZING  THE  KLEIN-GORDON  FORM  OF 
THE  STRATIFIED  A-LAYER  MODEL 

A  standard  second- order  central- difference  scheme  is  used  to  discretize  (V.19). 
On  a  spatial  and  temporal  grid  of  our  choosing,  we  let  be  the  FD  approximation 
of  7]i(x,y,t)  at  grid  point  (jp,  yq)  and  at  time  tn  in  layer  L*.  Thus  solving  explicitly 
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fox  tj?^1  yields: 


where  Q)t  =  \/g®i .  We  use  this  inter ioi  scheme  in  the  numerical  experiments  pre¬ 
sented  in  the  next  section.  The  Higdon  NRBC  is  applied  independently  to  each  layer 
and  takes  the  form: 

(V.21) 

The  development  of  its  discretization  follows  the  discussion  in  Section  IV.C.  Hard  wall 
or  Neumann  conditions  are  also  applied  independently  to  each  layer  and  discretized 
accordingly.  Since  the  discrete  forms  of  the  Higdon  NRBC,  Nuemann  condition,  and 
the  time  stepping  scheme  (V.20)  are  explicit,  the  whole  scheme  is  explicit. 


C.  NUMERICAL  EXAMPLE:  THE  STRATIFIED  Y- LAYER 
MODEL 

Consider  a  geophysical  process  that  occurs  in  a  semi-infinite  channel  (Fig¬ 
ure  9).  All  assumptions  and  simplifications  used  to  derive  the  TV-layer  model  apply. 

A  Cartesian  coordinate  system  (x, y)  is  introduced  such  that  the  channel  is  parallel 
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to  the  x- direction.  On  the  north  and  south  boundaries  Tjv  and  Ts  we  specify  the 


Neumann  condition: 

=  0  foi  i  =  1...7V. 

On  the  west  boundary  rH,-  we  prescribe  using  a  Dirichlet  condition,  i.e., 


(V.22) 


V t(°.  y.  *)  =  Viv.  (y,  0  foj  *  = 


(V.23) 


where  77, Vi  (1/,  f)  is  a  given  function  for  an  incoming  wave  or  disturbance.  At  x  — *•  00  the 
solution  is  bounded  and  does  not  include  any  incoming  waves.  The  initial  conditions 
are: 

J7*(x,  y,  0)  =  0,  =  0  for  i  =  l ...N.  (V.24) 

To  obtain  a  well-posed  problem  in  a  finite  domain  ft  we  impose  a  high- order  Higdon- 
NRBC  on  the  east  boundary  r^. 

We  now  apply  the  new  stratification  scheme  to  a  test  problem  using  the  semi- 

JY 

infinite  wave-guide.  A  channel  width  6  =  5  and  depth  =  .1  are  se- 

t=L 

lected  (note  that  to  satisfy  the  shallow  water  assumption,  it  must  be  true  that  depth 
d  ^horizontal  dimensions).  The  stratified  medium  is  modeled  with  six  layers.  The 
layer  thicknesses  from  top  to  bottom  are  =  {.01,  .01,  .01,  .01,  .01,  .05} .  The  density 
for  each  layer  is  given  by  p.  =  (1, 1.05,  1.1, 1.15, 1.2, 1.25} .  A  gravitational  parameter 
g  =  10  and  a  dispersion  parameter  /  =  .5  is  used. 

The  boundary  function  on  rH.-  is  stipulated  to  simulate  two  geophysical 
events.  A  surface  disturbance,  akin  to  the  wind  acting  on  the  ocean,  is  initiated  in 
Li  by  setting: 


.  ,v  f  .0005  cos  —  (y  —  2.5)  sinTrf  0  <  t  <  12.5  ,  ... 

7«Ay>*)  =  \  L5'y  '}  -  -  (V.25) 

\  0  otherwise, 


Note  that  the  maximum  amplitude  of  the  disturbance  is  small  relative  to  the  layer 
thickness  dl  =  .01  so  that  the  validity  of  the  model,  which  is  based  on  perturbation 
analysis,  is  not  violated.  A  second  disturbance,  simulating  seismic  activity  on  the 
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ocean  floor,  is  initiated  in  the  Lq  and  given  by: 


= 


if  |y-  2.5|  <  1.5  L  6  <  f  <  7.5 
otherwise  . 


(V.26) 


Ail  other  values  for  rjn  ,  are  zero.  The  simulation  is  run  for  15  time  units. 

The  problem  is  solved  for  three  different  scenarios.  First,  an  extended  domain 
V  is  constructed  using  a  15  x  5  rectangle  with  a  60  x  20  mesh  to  compute  a  reference 
solution  rjref.  Then  two  additional  solutions,  and  qmml,  are  computed  on  a 

truncated  domain  ft  in  which  an  artificial  boundary  B  imposed  at  x  =  5.  A  20  x  20 
mesh  is  used  on  the  resulting  5x5  square  so  that  the  mesh  for  ft  and  V  are  identical 
in  the  truncated  region,  For  ,  a  Higdon  NRBC  of  order  J  =  5  with  parameters 
Cj  =  {1,  1, 1, 1, 1}  is  used  for  each  layer.  For  qsa„2,  J  =  2  and  Cj  =  {1, 1}  is  used. 
Note  that  the  value  q  in  each  case  represents  the  total  perturbation  on  the  domain 
surface  and  is  a  combination  of  perturbation  effects  in  each  layer, 

The  reference  solution  q^f  is  then  juxtaposed  with  and  q _ ,  both  graph¬ 

ically  and  quantitatively.  The  respective  numerical  solutions  are  used  to  produce  filled 
contour  plots  that  provide  a  “fingerprint75  for  the  resulting  wave  action.  All  contours 
and  color  schemes  are  relative  to  qref.  These  solutions  are  then  used  to  obtain  error 


measurements  II 0 


which  are  calculated  by: 


=  hh\  N*N* 


(V.27) 


where  Nx  and  jYy  are  determined  by  grid  spacing,  In  both  cases,  the  error  at  the 
surface  and  for  each  layer  interface  is  plotted  versus  time.  Note  that  the  size  of 
V  precludes  spurious  reflections  from  polluting  q^f-  Therefore  ||e(t)||n  serves  as  a 
measure  of  spurious  reflection  at  B. 

At  t  =  7  (Figure  28),  the  surface  disturbance  in  has  populated  ft  and 
is  now  visible  in  I?.  In  addition,  the  bottom  disturbance  has  been  initiated  and  is 
propagating  in  ft,  Both  1  and  qsml  exhibit  wave  traces  similar  to  those  in  q^f . 
However,  the  error  measurement  for  q^^  is  an  order  of  magnitude  smaller  than  that 
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Figure  28.  HNRBO2DR-1S-6L-U0-V0-T07 


□f  T7„„, .  This  demonstrates  that  simply  increasing  the  Higdon  NRJ3C  order  J  reduces 
spurious  reflection.  A  description  of  the  shorthand  notation  used  in  the  figure  caption 
was  given  in  Section  IV. E. 

At  time  t  =  15  (Figure  29),  the  surface  disturbance,  which  ended  at  t  =  12.5, 
continues  to  propagate  in  and  V.  The  bottom  disturbance  has  successfully  passed 
through  rE  without  a  significant  increase  in  spurious  reflection.  The  wave  trace  for 
closely  resembles  that  of  Tfref,  however  deviations  in  are  now  visible.  This 
example  demonstrates  that  a  properly  constructed  Higdon  NRBC  can  be  used  to 
restrict  the  domain  of  the  N-layer  stratification  model  governed  by  the  linearized 
shallow  water  equation. 
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Figure  29.  HNRBO2DR-1S-6L-U0-V0-T15 
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D.  ADDITIONAL  OBSERVATIONS  FOR  THE 
STRATIFIED  V-LAYER  MODEL 

Using  the  same  problem  parameters  as  in  the  previous  example,  we  now  con¬ 
sider  the  behavior  of  the  perturbation  at  each  layer  as  predicted  by  the  model.  Filled 
contour  plots  are  constructed  to  measure  the  magnitude  of  perturbation  for  the  layers 
Ll  through  L§ .  All  contours  and  color  schemes  are  relative  to  the  perturbation  in  L$ . 
In  Figure  30,  we  see  that  the  action  for  both  the  surface  and  bottom  initiated  events 
has  affected  all  layers,  but  most  of  the  energy  transmitted  by  the  wave  “suiks”  to  the 
lower,  denser  layer  Lq.  This  remains  true  at  t  =  15  (Figure  31)  at  which  time  both 
events  have  passed  through  the  truncated  domain  12. 
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Figure  30.  HNRBC-2DR-1S-6L-U0-V0-T07  Layer  Perturbation  Comparison 


We  further  consider  model  predictions  for  the  behavior  of  the  perturbation 
at  each  layer  interface.  Recall  that  the  disturbance  at  a  layer  interface  is  the  sum 
of  the  disturbances  of  all  layers  below  the  interface.  Filled  contour  plots  are  now 
constructed  for  the  surface  and  the  five  layer  interfaces  using  data  obtained  from 
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Figure  31.  HNRBG-2DR-1S-6L-U0-V0-T15  Layer  Perturbation  Comparison 


AH  contours  and  color  schemes  are  relative  to  the  total  surface  perturbation. 
The  first  comparison  is  again  at  t  =  7  (Figure  32).  The  contour  plots  reveal  that  the 
magnitude  of  the  response  to  the  surface  event  is  damped  with  each  successive  layer 
interface  downward.  However,  the  wave  appears  to  maintain  its  character  at  each  layer 
interface  as  far  as  wave  speed  and  geometric  dispersion  is  concerned.  We  also  note 
that  the  bottom  layer  event  on  the  left  side  of  the  plot  has  immediately  propagated 
to  the  surface.  This  is  reasonable  since  the  fluid  is  incompressible.  However,  as 
the  wave  propagates  through  fi,  the  effect  on  the  lower  layer  interfaces  once  again 
dampens.  At  t  =  15,  the  damping  phenomenon  in  the  lower  layers  is  more  pronounced 
(Figure  33).  However,  there  is  an  exception  on  right  side  of  the  plot.  In  this  regime, 
the  perturbation  is  nearly  zero  and  affected  primarily  by  spurious  reflection  from  B. 
Differences  in  each  layer’s  reflection  may  have  caused  this  visible  anomaly. 

For  a  final  comparison,  we  consider  a  six-layer,  two-layer,  and  single-layer 
model.  For  this  we  set  up  the  parameters  for  two- layer  model  as  follows:  8i  = 
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Figure  32.  HNRBC-2DR-1S-6L-U0-V0-T07  Laver  Interface  Perturbation  Comparison 

{.05,  .05}  and  p.  =  {1.1, 1.25} .  Essentially  we  have  combined  the  upper  five  layers  of 
the  six- layer  model  and  used  the  average  density  for  the  first  layer  of  the  two- layer 
model.  The  second  layer  of  the  two- layer  model  is  the  same  as  the  bottom  layer 
of  the  six-layer  model.  The  parameters  for  the  single-layer  model  are  6  =  .1  and 
p  =  1.175  (note  that  from  (V.19)  we  conclude  that  density  is  not  a  parameter  in  the 
single-layer  model  and  therefore  the  value  for  p  is  irrelevant  ).  The  problem  is  again 
run  for  15  time  units,  however  the  domain  is  extended  to  x  =  15  for  each  model  and  a 
60  X  20  mesh  is  used.  This  eliminates  artificial  boundary  effects  and  comparisons  can 
be  made  without  concern  for  spurious  reflection.  All  other  parameters  are  the  same. 
The  results  of  the  three  models  are  presented  for  t  =  15  (Figure  34).  Comparing 
the  contours  of  the  models  reveal  that  the  surface  effect  is  reduced  as  the  number  of 
layers  increases.  Hence,  the  single-layer  representation  tends  to  over-predict  surface 
wave  action  if  the  medium  is  stratified. 

We  now  wish  to  apply  the  TV-layer  model  to  an  ocean  environment.  However, 
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Figure  33.  HNRBC-2DR-1S-6L.-U0-V0-T15  Layer  Interface  Perturbation  Comparison 


before  doing  so,  we  must  first  extend  the  application  of  the  Higdon  NRBC  to  two 
dimensions.  This  will  enable  us  to  model  a  “patch  of  ocean”  vice  “a  semi-infinite 
channel  with  hard  walls” . 
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Figure  34.  Model  Comparison:  HNRBCC2DR-1S-6L-U0-V0-T15  vs.  HNRBC-2DR- 
1S-2L-U0-V0-T15  vs,  HNRBO2DR-1S-1L-U0-VQ-T15 
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VI.  APPLYING  HIGDON  NRBC’S  TO  TWO 
OR  MORE  SIDES  OF  A  TWO-DIMENSIONAL 
RECTANGULAR  DOMAIN 


The  equations  and  concepts  are  now  in  place  to  construct  finite  domains  in 
which  Higdon  NRBC’s  are  imposed  at  two  or  more  boundaries.  In  this  chapter 
we  employ  the  use  of  several  examples  to  check  the  viability  and  the  stability  of 
schemes  that  utilize  Higdon  NRBC’s  on  two  and  four  sides  of  the  single-  and  multi¬ 
layer  models.  Special  attention  is  given  to  the  corners  where  two  Higdon  boundaries 
intersect. 


A.  EXAMPLE  ONE:  A  TWO-DIMENSIONAL  SINGLE¬ 
LAYER  SCHEME  WITH  HIGDON  NRBC’S  ON  TWO 
SIDES 

For  this  example  we  employ  a  quarter  plane  shown  in  Figure  35.  The  semi- 
y 

D 


Figure  35.  The  Semi-Infinite  Quarter- Plane 


infinite  domain  V  is  bounded  by  T s  and  rH/  and  is  represented  by  a  10  x  10  square 
with  a  40  x  40  mesh.  The  truncated  domain  is  a  5  x  5  square  with  a  20  x  20  mesh 
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ajid  bounded  by  Tat,  Ts,  r^,  and  T «. .  On  both  domains  Ax  =  .25,  Ay  =  .25,  and 
At  =  .025.  A  gr avitatdonal  parameter  g  =  10  and  a  dispersion  parameter  /  =  .5  is 
used.  A  single-layer  non- stratified  medium  with  a  thickness  of  d  —  0.1  is  considered. 
The  problem  is  run  for  15  time  steps. 

The  specification  for  the  Higdon  NRBC  on  Tjy  is  given  by: 


n(s+<*£)’-°- 


(VI. 1) 


This  equation  is  the  r/- direction  analog  of  (IV. 9)  which  was  used  to  describe  the 
Higdon  condition  for  r^.  A  discussion  equivalent  to  that  in  Section  III.C  applies  in 
discretizing  (VI. 1).  If  Aj  =  Ay  on  the  specified  grid,  and  the  same  Cj's  are  used  for 
each  boundary,  then  the  parameters  and  operators  used  to  discretize  the  rr  boundary 
condition  may  also  be  used  to  discretize  the  Tjv  boundary  condition.  To  complete 
problem  description,  a  Neumann  condition  is  imposed  at  T s  and  given  by: 


and  a  boundary  function  is  imposed  on  fV  and  is  given  by: 


H  Ar 


cos 


'nmK{y  -  1.875)  \ 


sm(^f) 


fT3=l 

0 


3.75  j 

where  the  parameters  A^,  and  n>m  are  selected  a  priori  as  follows: 


if  y  <  3.75, 
otherwise  , 


(VI  .2) 


Am  =  .001,  .002,  .001: 

=  1,  3,  1: 

u>m  =  .81,  1.37,  1.68, 


(VI. 3) 


The  restriction  y  <  3.75  is  discussed  in  the  next  section.  Homogeneous  initial  condi¬ 
tions  apply. 

A  solution  t -jrej  and  is  computed  for  V  and  respectively.  Filled  contour 
plots  are  generated  for  each  solution.  In  Figure  36,  these  solutions  are  reported  at 
t  =  3.  The  upper-right  subplot  displays  the  solution  for  on  V.  The  upper-left 
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Figure  36.  HNRBO2DR-2S-1L-U0-V0-T03:  Wave  Pulse  Passes  through  T* 


subplot  is  a  magnification  of  77  on  the  sub-domain  ft  and  represents  a  solution  on  ft 
as  if  no  boundaries  were  present  atJ  =  5andT/  =  5.  The  lower-left  subplot  reports 
the  solution  for  rjcasei  011  ^  where  Higdon  NRBC’s  are  imposed  at  and  A 
qualitative  comparison  can  be  made  between  the  two  left-side  subplots  to  determine 
the  effectiveness  of  the  NRBC  in  restricting  the  domain.  A  quantitative  measurement 
of  this  comparison  appears  in  the  lower-right  subplot,  which  reports  an  error  measure 
lk(0lla  (See  IV. 37).  Note  that  for  computational  purposes,  Higdon  NRBC’s  with 
Cj  =  (1, 1, 1, 1, 1, 1, 1}  are  imposed  at  y  =  10  and  r  =  10  on  P,  Any  reflection  from 
these  extended  boundaries  of  V  will  have  little  or  no  effect  on  the  domain  of  interest 
ft. 

In  Figure  36  at  t  =  3,  the  wave  pulse  has  reached  Tat.  A  visual  inspection 
reveals  no  discernible  spurious  reflection.  The  magnitude  of  the  error  measurement 
is  within  acceptable  norms  (i.e.  less  than  the  error  due  to  the  discretization). 

Atf  =  5  (Figure  37),  the  wave  pulse  has  reached  and  passes  without  visible 
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Figure  37.  HNRBa2DR-2S-lL-U(hV0-T05:  Wave  Pulse  Passes  through  VE 


reflection  or  undue  increase  in  ||e(f)||n.  Note  that  a  leading  edge  of  the  pulse  is  about 
to  reach  the  intersection  of  Tjv  and  At  t  =  15  (Figure  38),  the  wave  has  passed 
through  the  corner  unimpeded  with  no  adverse  effect  on  ||?(t)||a  and  the  two  left¬ 
side  subplots  are  similar  with  no  visible  differences.  The  measured  error  ||e(t)||n  has 
stabilized  at  approximately  2  x  10-6. 

The  corner,  however,  is  still  cause  for  concern,  because  two  different  ways  to 
approximate  the  Higdon  boundary  grid  points  are  possible,  In  the  first  method  we 
initially  evaluate  the  Te  grid  points  excluding  the  corner  point.  We  then  evaluate 
the  Tjv  grid  points  including  the  corner  point.  Hence  the  corner  point  is  evaluated 
based  on  boundary  values  calculated  for  Te.  This  method  was  used  in  generating 
Figures  36  through  38.  In  a  second  method,  the  grid  points  for  Tjv  (excluding  the 
corner  point)  are  fully  evaluated  first  followed  by  the  grid  points  of  Te  (including  the 
corner  point).  Now  the  corner  paint  values  are  obtained  from  T N  boundary  values. 
Regardless  of  the  method  used,  the  value  of  the  corner  point  should  be  the  same.  This 
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Figure  38.  HNREC-2DR-2S-1L-U0-V0-T15:  Wave  Pulse  Passes  through  Corner 


was  tested  and  the  results  are  compared  via  the  error  measurement  plot  in  Figure  39. 
The  figure  indicates  that  the  two  solutions  are  exactly  the  same.  Hence  no  special 
considerations  must  he  given  when  handling  corner  points  that  are  the  intersection 
of  the  Higdon  boundaries.  However,  one  must  be  careful  to  ensure  that  the  corner 
point  is  not  evaluated  twice  when  designing  algorithms  to  compute  these  boundary 
values. 
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Figure  39.  Corner  Point  Check:  Two  Approaches  to  Evaluating  HNRBC  Grid  Points 


B.  INSTABILITIES  IN  THE  2-D,  2-SIDED  HIGDON  NRBC 
SCHEME 


Figure  40.  HNRBO2DR-2S-1L-U0-V0-T03:  TN  Instability  with  Cj  =  {1, 1, 1, 1, 1} 


Despite  the  heartening  results  of  the  previous  example,  a  stability  problem 
does  emerge  along  the  entire  northern  NRBC  if  a  non- zero  boundary  function  rjw 
is  extended  to  Recall  the  constraint  y  <  3.75  in  (VI .4).  This  was  specified  in 
order  to  set  up  a  buffer  consisting  of  five  zero  valued  grid  points  between  the  boundary 
function  in  the  northern  Higdon  boundary.  If  this  is  not  done,  the  boundary  instability 
occurs  and  propagates  through  fi.  For  example,  define  the  boundary  function  at  T «-• 
as: 

Vm:  (y>  *)  =  Am  COS  sintyW)  (VI  .4) 

The  grid  point  at  y  =  5  is  immediately  adjacent  to  rty .  All  other  parameters  with 
regards  to  the  previous  example  remaining  unchanged.  The  problem  is  run  with  J  =  5 
and  Cj  =  (1,1,1, 1,1}.  In  Figure  40  at  t  =  3,  we  discern  turbulence  in  the  upper 
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left-hand  corner  of  ft  for  rjcate !•  In  Figure  41  atf  =  10,  we  see  that  the  instability  has 
propagated  through  ft. 


Figure  41.  HNRBO2DR-2S-1L-U0-V0-T10:  TN  Instability  with  Cj  =  {1, 1, 1, 1, 1} 


To  explain  the  unstable  behavior,  recall  that  the  5^-  order  Higdon  NRBC  on 
Tjv  is  constructed  with  up  to  5^-order  derivatives  with  respect  to  y  and  t.  Using 
the  backward- difference  method  to  construct  the  discretized  forms  of  these  derivatives 
requires  that  we  reach  back  as  many  as  five  grid  points  both  temporally  and  spatially. 
During  the  initial  time-stepping  sequence,  we  only  have  a  history  of  one  time  step. 
The  current  algorithm  uses  a  default  of  zero  if  the  discretization  calls  for  data  at  a 
time  that  is  less  than  zero.  This  is  a  poor  estimate  if  the  initial  condition  is  some 
non-zero  value  because  it  introduces  a  discontinuity  in  time.  The  problem  corrects 
itself  somewhat  if  a  Sommerfeld  condition  ( J  =  1  and  Cj  =  {1})  is  utilized  on  Tat 
as  in  Figure  42.  However,  our  goal  is  to  use  high-order  Higdon  conditions  to  increase 
the  accuracy  of  the  estimate.  Experimentation  showed  that  a  buffer  of  J  grid  points 
is  necessary  to  achieve  stability  for  a  J^-order  Higdon-NRBC.  The  exact  nature  of 
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the  instability  must  be  further  investigated.  However,  for  the  time  being  we  continue 
to  present  examples  while  maintaining  the  required  buffer  zone  to  achieve  stability. 
This  is  one  of  the  reasons  Givoli  and  N eta  introduced  auxiliary  variables  [Ref.  24]. 
See  also  Givoli  et  al  [Ref.  26], 


Figure  42.  HNRBC-2DR-2S-1L-U0-V0-T10:  Mitigation  of  Instability  using  Sornmer- 
feld  Condition  ( J  =  1  and  Cj  =  {!}) 
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C.  EXAMPLE  TWO:  A  TWO-DIMENSIONAL  MULTI¬ 
LAYER  STRATIFIED  SCHEME  WITH  HIGDON  NRBC’S 
ON  TWO  SIDES 


In  this  section  we  check  the  compatibility  of  the  two-sided  Higdon  model  with 
the  multi-layer  stratification  model  presented  in  Chapter  5.  This  example  will  again 
use  the  semi-infinite  quarter-plane  (Figure  35)  and  the  same  domain  specification  as 
used  in  Example  1  of  this  chapter.  The  stratified  medium  is  modeled  with  six  layers. 
The  layer  thicknesses  from  top  to  bottom  are  8{  =  (.01,  .01,  .01,  .01,  .01,  .05} .  The 
density  for  each  layer  is  given  by  p.  =  (1, 1.05, 1.1,  1.15, 1.2, 1.25}.  A  gravitational 
parameter  g  =  10  and  a  dispersion  parameter  /  =  .5  is  used. 

The  boundary  function  rjn.  on  IV  is  stipulated  to  simulate  two  geophysical 
events.  A  surface  disturbance  akin  to  the  wind  is  initiated  in  L\.  Using  the  function: 

=  51  An  cas  ^L875^  srn(wmt),  (VI. 5) 

TT3=L  \  / 


(VI. 6) 


we  let: 

if  t  <  12.5  L  y  <  3.75, 

0  otherwise 

where  and  are  given  by  (VI.  3).  A  second  disturbance,  simulating  seismic 

activity  on  the  ocean  floor,  is  initiated  in  the  Lq  using  the  function: 


,  .  .  /  x(y  —  1.875)  \  / 7r(t  —  6.5) ' 

/!(j,,i)  =  .001cosf - — - I  cos( - - - 


(VI  .7) 


and  letting: 

if5<f<8fcl<2/<2.75, 

.  (VL8) 

1 0  otherwise 

All  other  values  for  r/% v,  are  zero.  The  simulation  is  run  for  15  time  units  and  results 
are  reported  using  a  similar  presentation  as  that  used  in  Example  1. 

At  t  =  6  (Figure  43),  the  surface  effect  has  passed  through  Tjv  and  and 
||e(t)||n  at  the  surface  falls  within  acceptable  tolerances.  The  effects  of  the  bottom 
disturbance  are  readily  visible  near  the  west  boundary.  At  t  =  15  (Figure  44),  the 
bottom  disturbance  has  propagated  through  and  has  successfully  passed  through 
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Figure  43.  HNRBa2DR-2S-6L-U0-V0-T06  with  J  =  5  arid  Cj  =  {1, 1, 1, 1, 1}  fox 
each  Layer 


both  artificial  boundaries.  Some  differences  between  and  rjcasei  aze  visible,  but 
again,  the  error  measure  ||e(t)||n  is  less  than  the  discretization  error.  From  this 
example,  we  conclude  that  the  multi-layer  stratification  model  will  work  effectively 
when  Higdon  NRBC’s  are  placed  on  two-sides  of  the  domain 
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Figure  44. 
each  Layer 


HNRBC-2DR-2S-6L-U0-V0-T15  with  J  =  5  and  Cj  =  {1, 1, 1, 1, 1}  for 
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D.  EXAMPLE  FOUR:  A  TWO-DIMENSIONAL  MULTI¬ 
LAYER  STRATIFIED  SCHEME  WITH  HIGDON  NRBC’S 
ON  FOUR  SIDES 

In  this  example  we  consider  a  truncated  domain  where  Higdon  NRBC’s  are 
imposed  on  four  sides.  The  boundary  conditions  on  Tjv  and  Re  are  discretized  as 
before.  The  boundary  conditions  on  Ts  and  are  discretized  in  a  similar  manner 
(see  Section  III.C),  however  forward-  vice  backward- difference  equations  are  used 
in  approximating  the  spatial  derivatives.  As  before,  SI  is  a  5  x  5  square  with  a 
20  X  20  mesh.  Higdon  NRBC’s  are  located  at  x  =  0,  y  =  0,  x  =  5,  and  y  =  5  with 
Cj  =  {1, 1, 1, 1, 1}  at  each  boundary. 

The  extended  domain  V  is  the  infinite  plane.  It  is  represented  by  a  15  X  15 
square  with  a  60  x  60  mesh.  The  domain  of  interest  R  is  located  in  the  center  of 
V  at  5  <  x,y  <  10.  Higdon  boundaries  are  placed  at  x  =  0,  y  =  0,  j  =  15,  and 
=  15  with  Cj  =  (1, 1, 1, 1, 1, 1, 1}  at  each  boundary.  This  is  done  for  computational 
purposes  only  and  any  spurious  reflection  from  these  extended  boundaries  will  not 
significantly  pollute  the  domain  of  interest.  On  both  domains  Ax  =  Ay  =  .25  and  the 
time- step  At  =  .0125.  A  gravitation  parameter  of  g  =  10  and  a  dispersion  parameter 
of  /  =  .5  is  used.  The  layer  thicknesses  are  6{  =  {.01,  .01,  .01,  .01,  .01,  .05}  with 
respective  densities  given  by  p.  =  {1, 1.05, 1.1, 1.15, 1.2, 1.25}. 

Random  values  that  represent  physical  disturbances  in  are  introduced  at 
designated  times  at  specified  x-  and  ^/-ranges  in  specified  layers.  At  the  given  time, 
these  random  values  are  added  to  existing  values  of  rj  in  the  designated  layer.  The 
dispersive  mechanism  defined  by  (11.88)  is  than  allowed  to  act  on  the  event.  These 
random  events  are  sufficiently  complex  in  nature  and  “flex7’  the  Higdon  NRBC’s 
exposing  any  unacceptable  spurious  reflection. 

For  this  example  two  separate  disturbances  or  events  are  imposed.  Event  1  is 
a  surface  event  given  by: 

+  ( .0001  *  ran<f(— .  5,  .5)  if  1.5  <  x,  y  <  3.5, 

S^026  =  1  K  “  (VI. 9) 

L  0  otherwise 
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where  S^  025(x,y)  represents  a  disturbance  initiated  in  the  first  layer  Li  at  t  =  .025 
and  rand{— .5,  .5)  is  a  landom  number  on  the  interval  [—.5,  .5],  Event  2  is  a  bottom 
event  given  by: 

f. 00015  *  r<m<*(-.25, . 75)  if  1.5  <  x  <  2.25  L  1.5  <  y  <  3.5,  , 

Sj=5  =  i  K  '  “  “  _  (VI. 10) 

l0  otherwise. 

where  5^5(x,j/)  indicates  that  the  event  was  initiated  in  L$  at  t  =  5.  The  maximum 
amplitude  of  each  event  is  set  to  be  a  small  fraction  of  the  total  thickness  of  the  layer 
in  which  it  is  initiated,  thus  maintaining  the  integrity  of  the  perturbation  analysis 
that  was  utilized  to  derive  the  Klein-Gordon  equation.  A  buffer  of  at  least  5  zero- 
valued  grid  points  was  maintained  between  the  NRBC  and  the  event  acknowledging 
the  stability  considerations  discussed  in  Section  VI.  C.  All  events  must  be  shifted  5 
units  in  the  x-  and  y-directions  when  initiating  activity  on  V,  in  order  to  properly 
place  them  in  the  domain’s  center ,  The  random  seed  for  these  events  are  preserved 
for  use  in  later  examples. 
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Figure  45.  HNRBC-2DR-4S-6L-U0-V0-T01:  Surface  Disturbance  Initiated. 
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Figure  46.  HNRJ3C-2DR-4S-6L-U0-V0-T06:  Surface  Disturbance  has  Left  ft.  Bottom 
Disturbance  Initiated. 


A  trial  is  run  for  15  time  units.  At  t  =  1  (Figure  45),  the  surface  event 
has  been  initiated  and  the  resulting  waves  are  propagating  toward  the  four  sides  of 
ft.  The  error  measurement  ||e(t)||n  is  near  zero.  At  t  =  6  (Figure  46),  the  surface 
event  has  left  ft  with  no  significant  spurious  reflection  at  either  the  boundaries  or  the 
corners.  Additionally,  the  bottom  event  has  been  initiated  and  is  now  propagating 
outward  in  ft.  At  t  =  13  (Figure  47),  both  events  have  passed  through  the  artificial 
boundary  successfully.  An  outward  radiating  wave  front  can  be  seen  in  the  upper- 
right  subplot  of  the  extended  domain  V.  The  left-side  sub-plots  are  remarkably 
similar  indicating  that  the  boundary  conditions  are  suitable  for  the  layered  problem. 
At  t  =  15  (Figure  48),  the  wave  energy  resulting  from  the  disturbances  have  left 
ft.  The  effect  due  to  spurious  reflection  starts  to  become  significant  thus  increasing 
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Figure  47.  HNRBC-2DR-4S-6L-U0-V0-T13:  Ail  Events  have  Passed  through  HN- 
RBC’s, 
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Figure  48.  HNRBC-2DR-4S-6L-U0-V0-T15:  Spurious  Reflection  Predominates  in  f2. 
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VII.  LINEARIZED  SHALLOW  WATER 
EQUATIONS  WITH  NON-ZERO  ADVECTION 


In  the  this  chapter  the  linearized  shallow- water  model  is  further  developed  by 
lifting  the  zero-advectdon  r  equirement.  Rje call  the  linearized  form  of  the  shallow  water 
model  is  described  by  (11.77  and  11.79): 


du*  T1du*  T.du*  .... 


dv*  dv*  dv* 


St  +va~x  +  'ai  +  H°{aI  +  *7. 


-9 


dhB  Ago  dq\ 
dx  dx  dx  J 


J 


/  dkB  dH0  drA 

[dy  dy  dy) 


(VII. 1) 


0, 


where  U  and  V  are  constant  advection  terms  in  the  x-  and  ^-directions  respectively. 
As  per  discussion  in  Section  II. B. 8,  if  the  bottom  contour  hB  is  flat,  then  H0  is 
constant  .  This  allows  simplification  of  (VII .1)  to: 


cn 

1 

II 

II 

1 

‘cj 

(VII. 2) 

at  ox  ay  \  ox  ay  , 


=  0. 


Now  define  the  operator: 


D_  _  d_  d_  d_ 

Dt  dt^  dx  ^  dy' 


(VII. 3) 
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Substituting  this  into  (VII. 2)  yields: 


=  -9- 


:(0  +  K  +  /t ' 


=  -g- 


(VII  .4) 


D  .  rr  (dv.*  9v*\ 

+  +  air)  "  °’ 

where  fV  and  fU  are  constant  products  of  the  Cor iohs  parameter  and  the  advection 

components.  Note  that  this  is  analogous  to  the  zero- advection  form  of  the  shallow 

D  d 

water  equation  given  by  (11.78)  and  (11.80)  with  —  replacing  7-.  Equation  (VII  .4) 
is  combined  into  a  single  form  using  the  same  steps  delineated  in  Section  II, B. 9,  The 


result  is: 


D  (  D2  . 


Dt  V  Dt2 


fa)-CoV^  +  A  =0, 


(VII. 5) 


where  C0  =  y/gH0.  This  is  the  Klein- Gordon  form  of  the  linearized  shallow  water 
equation  with  non- zero  advection.  This  can  be  rewritten  as: 


-  clv2v  +  f2v  = 


(VII. 6) 


such  that: 


d  «  ds  as  ,.as  _ 


We  shall  consider  (VII. 6)  in  its  homogeneous  form: 

|£<>j>-c8v*,+a-o. 


(VII. 7) 


(VII. 8) 


Applying  the  operator  listed  in  (VII. 3)  twice  yields: 

D2  d2  ,  d2  ,  d2  d2  d2 

Dt2  dt2  dx2  dy 2  dxdt  dycft  dxdy 


(VII. 9) 
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Substituting  this  into  (VII. 8)  yields: 

(VII.10) 

This  is  an  expanded  Klein-Gordon  equivalent  for  the  linearized  shallow  water  equa¬ 
tions  with  non-zero  advection  terms  U  and  V .  It  applies  to  a  single-layer  model,  but 
will  later  be  extended  to  the  TV- layer  stratified  model. 


A.  DISCRETIZING  THE  LINEARIZED  SWE  WITH  CON¬ 
STANT  NON-ZERO  ADVECTION  TERMS 

The  following  central-difference  approximations  are  used  to  discretize  (VII.10): 

n +  o(Ax2) 


(VII.ll) 


TJxv~  4AxA  y 

Analogous  forms  are  used  for  approximating  rfa,  7fc.t  and  t/^.  Substituting  these 
into  (VII.10)  yields: 

^ +^L)  + (u 


Vp+l^+l  Vp+l.Q—l  Vp—l^+l  Vp-lj-l 


-|-  0(Ax2,  Ay2). 


Av2  ^p4+1  ^p-*?-1) 


(VII.12) 


l  ^  / ._rc  + 1  „n— 1  „rc+l  I  „n — L  v 

9 AvAf  vJ,p-<J+l  ^P4+1  ^p-?-1-  ^p^-i) 


£/V  , 


+  2AxAy^+l*+l  ~  ~  V^+l  +  =  °- 
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This  is  the  discrete  form  of  the  linear  shallow  water  equation  with  non-zero  advection 
terms.  Terms  that  contain  an  n  -f  1  superscript  are  unknown  rj  values  which  must  be 


solved.  Making  the  following  substitutions: 


1 

rj  (»’  -  Co) 

r  (v2-c2) 

At2’ 

Ax2  ’ 

Ay2 

V 

E_  V 

F-  UV 

2AxAt  ’ 

2  Ay  At  ’ 

2AxA  y’ 

and  moving  known  terms  to  the  right  side  of  (VII. 12)  yields: 


(VII. 13) 


ArC1  4  Dr%tl<  -  4  fig#,  -  fig#,  = 


(2-4  4  IB  4  2C  -  #),#  -  A£#  -  B  (>#.,„  4  <£_,«) 

(VII. 14) 

_  C  (jJpjj  +  L  +  ^  (^Pf  14  —  Vp-tg)  4-  E  —  Vpj-l) 


~  ^  (^+l-fl+L  —  Vp+ L.$— 1  Vp-  L$+l  Vp-ljj-l)  ' 

Since  there  are  five  unknowns,  (VI 1. 14)  must  be  solved  implicitly.  A  system  of  AT 
equations  with  N  unknowns  is  set  up  where  N  is  the  total  number  of  grid  points.  To 
make  the  indices  compatible  with  those  used  in  the  discretized  form  of  the  Higdon 
NRBC  (see  Sections  IV. C  and  IV.D)  we  replace  n  with  n  —  1.  This  yields: 


(VII. 15) 


The  Higdon  NRBC  (HI  .4)  is  now  applied  to  the  boundaries  of  the  problem.  The 
unknowns  in  the  discretized  form  of  Higdon  NRBC  (IV.  19)  are  those  terms  whose 
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operator  Pm  does  not  contain  a  time  shift  St  ,  For  example,  consider  the  discretized 
form  of  a  Higdon  NRBC  with  order  J  =  2  on  the  eastern  boundary  (IV.  18): 


[ai<r2  -(-  aid^Sf  -|-  aie2Sx  -)-  d]_a2St  -|- 

4-  die2SfS~  -1-  ei a2S~  -f  eid2SfSx  4-  e-i^S”2]^  =  0. 


(VII. 16) 


Shifting  the  known  quantities,  e.g,  ail  terms  with  time  shift  operators,  to  the  right 
side  of  (VII. 16)  yields: 


aid 2  4-  e2Sx  -(-  eia2Sx  4-  el e2Sx  2]rfeq  = 


-  [a]d2St  4-  dia2St  -|-  did2Si  2  -1-  die2St  Sx  -J-  eid2St  Sx ]rfeq, 
or  in  alternate  form: 


(VII. 17) 


<*L  *7nE.q  +  (al*2  +  +  eleiVE-74  = 


(VII. 18) 


-  (aLrf2  -|-  d^rfej  -  (die2  +  ei d2)rQ:^jq  -  d^rfe2. 


In  general,  the  number  of  unknowns  in  a  J^-order  Higdon  NRBC  using  first-order 
difference  approximations  is  J  -(- 1.  The  number  of  unknowns  in  a  J^-order  Higdon 
NRBC  using  second-order  difference  approximations  is  2J  -|-  1.  Equation  (VII. 15) 
together  with  an  analogous  form  of  (VII. 18),  adjusted  for  the  Higdon  order  J,  is  used 
to  set  up  the  system  of  equations  that  will  be  solved  to  “tame-step7’  the  system. 


B.  NUMERICAL  EXAMPLE:  TWO-DIMENSIONAL  SINGLE¬ 
LAYER  SCHEME  WITH  HIGDON  NRBC’S  ON  FOUR 
SIDES  WITH  NON-ZERO  ADVECTION 

In  this  example,  the  truncated  domain  12  as  described  in  Section  VI.  D  with 
Higdon  NRBC’s  on  four  sides  is  used.  As  before,  the  extended  domain  V  is  an  infinite 
plane  represented  by  a  15  X  15  square  with  a  60  X  60  mesh.  i2  is  located  in  the  center 
of  V  at  5  <  x,y  <  10.  Higdon  boundaries  are  also  imposed  on  V  for  computational 
purposes.  Spurious  reflection  from  these  boundaries  should  not  significantly  pollute 
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f2.  On  "both  domains  Ax  =  Ay  =  .25  and  At  -  .0125.  A  gravitation  parameter  of 
g  =  10,  dispersion  parameter  of  /  =  .5,  and  a  single  layer  of  thickness  ©  =  .1  with 
density  p  =  1  is  used.  Advection  constants  of  U  =  .5  and  V  =  —.25  are  -utilized. 

Physical  disturbances  in  ft  are  initiated  via  two  separate  events.  Event  1  is 
given  by: 


->f=.025 


(VII. 19) 


.0001  *  rand(—  .5,  .5)  if  1.5  <  x,  y  <  3.5, 

^  0  otherwise 

where  Sf=  025(j,?/)  represents  a  disturbance  initiated  at  f  =  .025  and  rand(— .5,  .5)  is 

a  random  number  on  the  interval  [—.5,  .5],  Event  2  is  given  by: 


St=t  = 


.000015  *  rand(— .25,  .75) 


if  1.5  <  x<  2.25  L1.5<y<  3.5, 


1 0  otherwise. 

(VII. 20) 

where  St=6(x,  y)  represents  a  disturbance  initiated  at  t  =  5  and  rand(— .25,  .75)  is 
a  random  number  on  the  interval  [—.25,  .75].  Note  that  these  events  are  analogous 
to  the  events  described  in  Section  VI. D.  The  random  seed  was  preserved  so  that  the 
same  values  would  be  generated  for  each  event.  A  buffer  of  at  least  5  zero- valued  grid 
points  was  maintained  between  the  NRBC  and  each  event  for  stability  purposes.  The 
events  are  shifted  5  units  in  the  positive  x-  and  y-  directions  on  I?  in  order  to  properly 
place  them  in  the  domain’s  center. 

Before  running  an  example,  consideration  was  given  to  the  selection  of  Cj' s. 
Several  experiments  were  conducted  with  results  reported  in  Figure  49.  Initially,  a 


Higdon  NRBC  with  order  J  =  5  and  Cj  =  {Co,  Co,  Co,  Co,  Co}  where  Co  =  y/g®  =  1 
was  considered.  This  was  compared  to  a  case  where  the  C/ s  are  corrected  for  advec¬ 
tion,  The  predominate  speed  of  the  gravity  wave  is  C0.  This  is  affected  somewhat 
the  dispersion  and  wave  height.  However,  with  the  inclusion  of  advection,  the  pre¬ 
dominate  wave  speed  with  respect  to  each  boundary  is  affected  more  significantly. 
Therefore  the  C^’s  on  each  boundary  are  adjusted.  These  adjustments  were  made  as 
follows: 

C?*  =Cj+U,  C*ezt  =  Cj-U, 

*rfll  =  Cj  +  V,  Cf 1 ^  =  Cj-  V. 


(VII. 21) 
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For  this  example,  the  adjusted  Cj’s  are: 


=  {1.5, 1.5, 1.5, 1.5, 1.5},  CJ Mrfh  =  {.75,  .75,  .75,  .75,  .75}, 

CJ"*  =  {-5,  -5,  .5,  .5,  .5},  =  {1.25, 1.25, 1.25, 1.25, 1.25}. 

The  results  of  both  runs  show  a  significant  decrease  in  ||e(t)||n  to  about  10-3  att  =10. 
This  error,  however,  can  be  reduced  further. 


Figure  49.  Plot  A:  J  =  5  with  Cj  =  {1, 1, 1, 1, 1}  adjusted  for  advection  compared 
to  J  =  5  with  Cj  =  {1, 1, 1, 1, 1}  unadjusted  for  advection.  Plot  B:  J  =  3  with 
Cj  =  {1,1,  1}  compared  to  J  =  5  with  Cj  =  {1,1,1, 1,1},  both  cases  adjusted  for 
advection.  Plot  C:  J  =  3  with  Cj  =  {.8,  .9, 1}  compared  to  J  =  3  with  Cj  =  {1, 1, 1}, 
both  cases  adjusted  for  advection.  Plot  D:  Corner  check  for  J  =  3  with  Cj  =  {.8,  .9, 1} 
adjusted  for  advection 


Recall  from  Section  VI. C  that  a  buffer  of  J  zero- valued  grid  points  was  neces¬ 
sary  to  achieve  stability  for  a  J^-order  Higdon  NRBC.  When  advection  is  incorpo¬ 
rated  into  the  problem,  this  buffer  zone  moves  horizontally  toward  at  least  one  of  the 
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boundaries.  Therefore  the  buffer  is  compressed  with  respect  to  the  boundary  toward 
which  it  is  moving.  In  order  to  maintain  stability  we  must  either  increase  the  size  of 
the  buffer  zone,  or  reduce  the  order  J.  In  plot  B  (top  right)  of  Figure  49  a  5th-order 
Higdon  NRBC  is  compared  to  a  3r<1-order  Higdon  NRBC.  In  both  cases  the  C/s  are 
adjusted  for  advection.  In  this  example,  ||e(t)||n  is  reduced  by  an  order  of  magnitude 
to  about  10— 'L. 

One  further  adjustment  is  possible  to  reduce  ||  e(t ) || n .  Geometric  dispersion 
is  another  factor  in  the  boundaries  response  to  an  impinging  wave.  A  wave  striking 
normal  to  the  boundary  will  generally  have  a  wave  speed  that  is  approximately  C0. 
In  all  other  cases,  the  wave  speed  is  less  than  Co.  An  example  was  set  up  for  J  =  3 
in  which  Cj  =  {.3,  .9, 1}  with  the  reduced  values  taldng  into  account  the  geometric 
dispersion.  Adjusted  for  advection,  the  C/s  used  for  the  problem  are: 


C?“f  =  {1.3, 1.4, 1.5},  C?0™  =  {.55,  .65,  .75}, 
C*est  =  {.3,  .4,  .5},  Cf*h  =  {1.05, 1.15, 1.25}. 


(VII. 22) 


In  Plot  C  (bottom  left)  of  Figure  49  an  additional  reduction  in  ||e(i)||n  is  evident. 
Further  analysis  is  necessary  to  determine  how  to  best  adjust  Cj  values  for  geometric 
dispersion. 

The  question  of  the  corner  points  of  is  again  salient  in  the  advection  case, 
because  the  values  for  Cj  on  each  boundary  are  now  different.  Recall  that  there  are 
two  ways  to  approximate  the  boundary  values  when  numerically  solving  the  problem. 
Both  approaches  are  tested  here,  In  the  first  run  the  j-boundaries  were  computed 
first  (including  the  corner  points)  and  the  ^-boundaries  computed  next  (excluding 
the  corner  points).  In  a  second  experiment  the  procedure  was  reversed  and  corner 
points  were  included  in  the  y- boundaries.  Plot  D  (bottom  right)  of  Figure  49  reveals 
that  the  solutions  are  identical.  Hence,  as  concluded  earlier,  no  special  handling  at 
the  corner  points  is  necessary. 

With  these  results  in  mind,  Higdon  NRBC’s  of  order  J  =  3  with  Cj  =  {.3,  .9, 1} 
are  used,  With  U  =  .5  and  V  =  —.25,  the  adjusted  Cj’s  are  those  listed  in  (VII. 22). 
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A  trial  is  run  fox  10  time  -units.  At  t  =  1  (Figure  50),  event  1  has  been  propagating 
outward  in  ft  for  approximately  1  time  unit.  The  effect  of  advectdon  is  apparent 
as  the  propagation  of  the  gravity  wave  is  tending  toward  the  southeast  (i.e.  in  the 
<  .5,  —.25  >  direction).  The  leading  edge  of  the  wave  has  passed  through  the 
but  the  error  measurement  is  still  very  small.  At  t  =  2  (Figure  51),  event  1  has 
crossed  and  Later,  at  t  =  3  (Figure  52),  event  1  has  crossed  Tjv  and  T H, .  At 
t  =  5  (Figure  51),  most  of  event  1  has  left  ft.  We  note  some  spurious  activity  on  the 
western  boundary. 

At  t  =  6  (Figure  54),  the  waves  generated  by  event  2  are  approaching  Fe  and 
r w  ■  Event  1  has  passed  through  all  four  boundaries  relatively  unperturbed.  The  plot 
of  V  reveals  that  the  wave  front  continues  to  tend  toward  the  southern  and  eastern 
portion  of  the  extended  domain. 
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Figure  50.  HNRBC-2DR-4S-1L-U.5-Vm.25-T01:  Event  1  Initiated. 


At  t  =  10  (Figure  55),  the  second  event  has  passed  through  the  boundary.  The 
wave  propagation  pattern  continues  to  “drift”  in  the  direction  of  advection  as  revealed 
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Figure  51.  HNRBC-2DR-4S-1L-U.5-Vm.25-T02:  Event  1  Crosses  Ts  and  r£. 


by  the  upper-right  plot  of  I?.  Close  inspection  of  the  contours  reveal  spreading  where 
the  gravity  wave  is  traveling  in  the  direction  of  advection  and  compression  where  the 
gravity  wave  is  traveling  against  the  direction  of  advection.  In  the  latter  case,  this 
indicates  a  steeper  wave  front.  Since  the  gravity  wave  is  omni-directional,  this  effect 
varies  throughout  the  plot.  In  the  noise  of  spurious  reflection  is  now  visible. 

This  experiment  was  repeated  for  two  other  sets  of  values  for  U  and  V .  In 
the  first  variation  (Figure  56)  the  magnitude  of  the  advection  constants  were  lowered 
to  U  =  .4  and  V  =  —.15.  As  expected,  there  is  a  decreased  tendency  toward  the 
southeast,  Also  notable  is  a  reduction  in  the  error  measurement,  In  the  second 
variation  (Figure  57),  the  magnitude  of  the  advection  constants  was  increased  to  U  = 
.6  and  V  =  —.35.  The  tendency  to  the  southeast,  as  well  as  the  error  measurement, 
has  increased.  These  results  indicate  that  the  model  is  behaving  as  expected  with 
regards  to  the  rate  and  direction  of  advection.  However,  as  the  magnitude  of  the 
advection  constants  is  increased,  the  measured  error  will  also  increase.  In  the  current 
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Figure  52.  HNRBC-2DR-4S-1L-U.5-Vm.25-T03:  Event  1  Crosses  FN  and  Fw 


example,  the  magnitude  of  the  advect ion  is  4  to  7  times  greater  than  the  magnitude 
of  the  depth.  In  a  real  world  problem,  where  the  open  ocean  is  the  medium  of 
propagation,  advection  constants  are  expected  to  be  significantly  smaller. 
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Figure  53.  HNRBC-2DR-4S-1L-U.5-Vm.25-T05:  Event  1  Leaves  ft  with  Visible  Spu¬ 
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Figure  54.  HNRBO2DR-4S-lL-U.5-Vm.25-T06:  Event  2  Initiated. 
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Figure  55.  HNREC-2DR-4S-1L-U.5-Vm.25-T10:  The  Noise  of  Spurious  Reflection 
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Figure  57.  HNRBC-2D-4S-1L-U.6-Vm.35-T10:  End  of  Run 
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C.  LINEARIZED  SWE  WITH  CONSTANT  NON-ZERO 
ADVECTION  TERMS  EXTENDED  TO  THE  Y- LAYER 
MODEL 

Recall  the  linearized  form  of  the  shallow  water  equation  with  non-zero  advec- 
tdon  terms  for  a  AT-laver  stratified  model  (V.10): 


^  I  TJ  |  V  &***  ffV  |  rM-  n  ^  I  Y  ^ 

aT  +  f7,  aT  +  m  +t()  “  ~s  i£*  &  +  £  JIl  • 


w- +  ■ +• ■ + m -  -  (SfM + s  t )  ■ 


IT  +  v<  3*  +  V|  Sy  +  e‘  Iftt  +  Sy)-0' 


where  U{  and  V \  are  the  x-  and  y-  components  of  advection  in  the  z^-layer .  The  same 
techniques  as  those  described  in  Section  II. B. 8  are  used  to  generate  a  single  equation 
for  t]i,  the  perturbation  of  the  layer.  As  an  intermediate  step  a  result  that  is 
analogous  to  (V.19)  and  (VII. 6)  is  generated: 


D: 


Dt 2 


't-L 


(vt )  -  g^tV2 1 Y1  ~Vj  +  H  ^ )  +  f2v<  =  si(z>  y>  *)> 

L  j=l  P*  )=i 


N 


(VII. 23) 


where  5,(x,  y,  t)  is  the  source  function  for  the  z1^- layer  of  an  AT-laver  model.  The 
^  is  defined  as: 


operator 


LDiJt 


St(x,  y,  t )  must  satisfy: 

D 


[Dt 


d  rr  9  ,-9 

&+0,Yv,V 


,c,  ...  as,  as,  as,  _ 

(>St(x,  y,  f ))  ^  Uj  T  Vj  0. 


(VII. 24) 


(VII. 25) 
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Considering  the  homogeneous  form  for  (VII. 23)  and  expanding  yields: 


(VII. 26) 


This  is  the  linearized  shallow  water  equation  for  an  AT-layer  model  with  non-zero 
advection  terms. 
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D.  DISCRETIZING  THE  LINEARIZED  SWE  N- LAYER 
STRATIFIED  MODEL  WITH  CONSTANT  NON-ZERO 
ADVECTION  TERMS 

As  befoie,  -using-  the  central- difference  approximations  in  (VII. 11)  to  discretize 
(VII.  26)  yields: 

7^ (Vt.w  ~  4-  Vt.m  )  + 


+ 


+ 

+ 


+ 


(7?i.p+l.<?  2r7t.p j  4  ^.p-l.<j)  "I"  ^^2  “^-PJ  "I"  ^-P4-l) 


_ /„«+!■  —  r?n-1  —  r7n+1  4-  r?n_1  ^ 

?AlAt  ['k-P+l4  ^i-p+l  4  i 


/  • n  +  1  „tj— 1 L  tj+1  .  Ti — L  X 

■.Vt.pj?+1  ^/s.pjj+L  *7t.p.a— L  '  Vi.pj}—  l) 


2AyAt 


?AxAy^i'p+lJl+l  ^-P-U+i  +  ^.p— u-l) 


t-L 


Z  7  (iffp+L*  -  2^  4  i5p-i„) 

j=l  rt 


(VII. 27) 


-  (  ^| )  Z  (tf-w+i  -  2t£p?  +  ’Sm-i) 


=  0. 


This  is  the  discretized  form  of  the  IV-layer  linearized  shallow  water  equation  with 
constant,  non- zero  advection  for  the  i^1- layer.  All  terms  that  contain  the  n  T  1 
superscript  are  unknown  values  of  rji  which  must  be  solved,  Making  the  following 
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substitutions : 


(VII. 28) 


moving  known  teims  to  the  light  side  of  (VII.27),  and  shifting  the  time  index  yields: 


+  DtVi.p+l.q  ~  .p-l.q  pj}  +  l~  ^Vi.p^-l  ~ 


(2A  +  2 B,  +  2 C,  - 


-  s.  ('fp+u  +  -  c.  +  »?”p.ti) 

+  D,  (^-+!u  -  +  E<  (>tf£+1  -  jjS3_l) 

IT  /  T3  —  1  T3 —  1  J!-l  i  fl— 1  ^ 

—  r«  Vft.pfljJ+l  _  '/t.pfl^-l  _  Vi.p-ljtj+l  "t  ’li.p-lj-l) 


+  G| 


+  G< 


+  « 


+  i/t 


t-i 


E  7  (^.p+i.9  -  2^.ml + u) 

j=i  ^ 


E  (Clu  -  2,?"«  +  >SA«) 


U=t 


t-1 


E£(*J«-aSS  +  >Ci-0 

j=i 


E  (^.p^J+L  2t7j-P9  “I"  ^-P4-l) 

j=i 


(VII.29) 


As  in  the  single-layer  advection  case,  (VII.29)  must  be  solved  implicitly  foi  each  layer 
L{.  The  system  of  equations  is  completed  on  the  boundaries  using  the  discretized 
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Higdon  'boundary  equations  as  discussed  in  section  VII.B. 


E.  NUMERICAL  EXAMPLE:  A  TWO-LAYER  SCHEME 
INCORPORATING  ADVECTION 

In  this  example  the  domains  ft  and  V  as  described  in  Section  VII.B  with  identi¬ 
cal  positioning  of  the  Higdon  NRBC’s  are  utilised.  The  following  problem  parameters 
are  used: 

Ax  =  Ay  =  .25,  At  =  .1, 

9  =10-  /  =  -5, 

©i  =  {.03, 07},  ^  =  {1, 1.05}, 

Ui  =  {.025,  .025},  V]  =  {-.025,  -.025}, 

J=  5,  C,  =  {.6,.7,.8,.9,1}. 

Correcting  the  C}’ s  for  advection  yields: 


(VII. 30) 


Ceast  =  {.625)  >725j  .825,  .925, 1.025},  C”orth  =  {.575,  .675,  .775,  .875,  .975}, 

CfV£t  =  {.575,  .675,  .775,  .875,  .975},  C =  {.625,  .725,  .825,  .925, 1.025}. 

A  single  physical  disturbance  is  initiated  in  ft  and  is  given  by: 

L  f  .000001  *rand(—. 5,  .5)  if2<r,i/<3, 

S-Li  =i  . 

1 0  otherwise  , 

where  y)  represents  a  disturbance  initiated  in  Ll  at  t  =  .1  and  rand{— .5,  .5) 

is  a  random  number  on  the  interval  [—.5,  .5],  The  example  is  run  for  five  time  steps. 

At  t  =  1  (Figure  58),  the  disturbance  has  been  underway  for  approximately 

one  second.  Minimal  spurious  reflection  occurs  at  the  boundaries.  In  the  lower-right 

plot  two  additional  measurements  are  noted.  The  first,  “Max  Ref  Surf75  is 

measured  over  the  entire  run.  The  next,  “Max  Hell  Ratio75  is  given  by: 


Max  | |e| |Ratio  =  (VII. 31) 

|  *(  | max 

Since  both  are  maximums  extracted  from  the  data  generated  over  the  entire  run, 
they  will  not  change  with  time.  Both  are  used  in  the  next  section  to  compare  Higdon 
boundary  effectiveness. 
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At  t  =  5  (Figure  59),  most  of  the  wave  action  has  left  ft.  The  residual  action  in 
the  truncated  domains  are,  for  the  most  part,  similar.  There  is,  however,  some  visible 
difference  near  the  boundaries  resultant  from  spurious  reflection.  The  lower-right  plot 
reports: 

Max  1 1  e 1 1 FLati o  =  1.08%. 

That  is  to  say,  the  maximum  error  norm  ||e||  at  t  =  5  was  1.08%  of  the  As  we 

shall  see  in  the  next  section,  this  is  a  favorable  measurement  for  the  Higdon  NRBC. 
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Figure  58.  HNRBC-2DR-4S-2L-Up025-Vmp025-T01:  Disturbance  Initiated. 
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Figure  59.  HNRBC-2DR-4S-2L-Up025-Vmp025-T05:  End  of  Run,  Some  Noise  at 
Boundaries  of  Bottom  Left  Plot. 


119 


F.  VARYING  PARAMETERS  FOR  THE  TWO-LAYER 
SCHEME  WITH  ADVECTION 

Thus  far  in  this  dissertation,  the  truncated  domain  ft  has  been  numerically 
compared  to  its  infinite  analog,  the  extended  domain  V,  Unfortunately,  most  of  the 
computational  assets  are  invested  in  generating  a  solution  for  this  “pseudo-infinite7’ 
domain.  In  the  four-sided  Higdon  NRBC’s,  the  extended  domain  V  contained  nine 
times  as  many  grid  paints  as  $7.  Addition  of  advection  considerations  required  em¬ 
ployment  of  implicit  methods.  In  these  cases,  computer  storage  became  a  key  factor. 
For  example,  V  for  a  four-sided  problem  with  AxA y  =  .25  required  a  612  x  612  ma¬ 
trix  to  generate  a  solution  whereas  ft  required  a  much  smaller  matrix  of  dimensions 
212  x  212.  Obviously  sparse  matrix  procedures  would  alleviate  this  requirement,  but 
the  point  illustrates  that  maintaining  an  extended  domain  becomes  unwieldilv  if  the 
domain  of  interest  is  enlarged  or  the  grid  is  refined. 

In  the  following  subsections,  parameters  of  a  two- layer  advection  problem  are 
varied  and  the  resulting  measurement  of  “Max  Hel  |  Ratio”  (VII.31)  is  compared.  From 
these  comparisons,  conclusions  are  listed  that  predict  parameter  limits  that  destabilize 
the  Higdon  NRBC.  Knowledge  of  this  behavior  will  allow  us  to  drop  extended  domain 
comparisons  and  focus  attention  on  the  domain  of  interest. 

In  all  trials,  the  domain  ft  is  a  5  x  5  square.  The  grid  is  varied  as  specified  in 
each  set  of  trials.  The  following  event  is  used  to  initiate  action  in  all  trials: 

f  .000001  *  randi— .5,  .5)  if  2  <  x,  y  <  3, 

S£At  =  l  m  (VII. 32) 

1 0  otherwise  , 

where  is  an  event  generated  in  the  t^-layer  at  t  =  At  (i.e.  the  first  time 

increment).  Each  trial  is  run  for  5  time  units  and  all  C^’s  are  adjusted  for  advection. 
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1.  Varying  Cj ’s,  Ax,  and  Ay 

In  this  set  of  trials  the  Higdon  coefficients  Cj  were  varied.  The  following 
parameters  were  fixed: 

Hi  =  {.03,  .07},  pi  =  {1.,  1.05}, 

Vt  =  {.025,  .025},  Vt  =  {-.025,  -.025},  (VII.33) 

At  =  .1  Ax  =  Ay  =  .5 

The  value  for  1^1^^  was  determined  to  be  1.92  x  10-e  and  was  unaffected  by  the 
choice  of  C/s.  Rjesults  for  Ax  =  Ay  =  .5  are  posted  in  Table  IV  and  were  compared 

1 \&\  \max 

using  ‘ — ,  a  measure  described  in  Section  VII.E.  This  set  of  numerical  trials  was 

1 7/ 1  max 


Table  IV.  Varying  C^-s  with  Ax  =  Ay  =  .5 


||g||n.<rr 
|T?I  max 

1HI"" 

lt?l  nux 

{■1,6,1} 

3.47% 

{6, .5,1} 

3.83% 

{.1,  .5,1} 

3.53% 

{•7,  .8,9} 

3.97% 

{.6, .8,1} 

3.57% 

{.8,  .9,1} 

4.13% 

{■1,75,1} 

3.58% 

{111} 

4.21% 

{.025,  .6, .975} 

3.58% 

{1,1} 

4.23% 

{.4,7,1} 

3.61% 

{.1,1,1} 

5.99% 

{.1,.5,.9} 

3.63% 

{1,1,1,!} 

92.2% 

{-5,5,5} 

3.76% 

repeated  for  Ax  =  Ay  =  .25.  In  this  case  1^1^^  was  determined  to  be  1.41  X  10“^, 
a  change  that  reflects  the  nature  of  the  random  event  (VII. 32)  that  was  used  to 
generate  the  disturbance.  The  results  of  these  trials  are  reported  in  Table  V.  Several 
conclusions  are  drawn  from  these: 

•  Refining  the  grid  results  in  a  more  efficient  boundary  conditions  (e.g,  less 
spurious  reflections). 

•  Higher  order  Higdon  NRBC’s  tended  to  generate  more  effective  boundary  con¬ 
ditions. 

•  A  buffer  zone  as  discussed  in  Section  VI. B  and  VII, B  must  be  considered. 
This  limits  the  Higdon  order  J  that  can  be  used,  For  Ax  =  Ay  =  .5,  the  trial 
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Table  V.  Varying  Cj-s  with  Ax  =  Ay  =  .25 


9 

3 

{.6, .7, .8, .9,1} 

1.08% 

{ .  1,  .75,1} 

2,14% 

(.9, .925, .95, .975,1} 

1.13% 

(.5,  .5,. 5} 

2.15% 

{i,i, i.i.i} 

1.16% 

{1,  6,1} 

2.24% 

{.1,.5,.6„7,1} 

1.19% 

lilGEjSSi 

2.37% 

P,. 8,9,1} 

1.21% 

{  •1,  5,  9} 

2.37% 

{1,1,1,!} 

1.24% 

{1,1} 

2.37% 

{.8, .9,1} 

1.57% 

(.025,  .6, .975} 

2.43% 

{1,1,1} 

1.58% 

(0,  .5,1} 

2.64% 

{.7, ,8, .9} 

1.60% 

{.i,.i,i} 

3.71% 

{,6, ,8,1} 

1.63% 

liHMIM 

393% 

{.4, .7,1} 

1.77% 

became  unstable  for  J  =  4.  For  Ax  =  Ay  =  .25,  instability  occur ed  at  J  =  6. 
In  the  latter  case,  higher  order  Higdon  NRBC’s  could  be  used  increasing  the 
effectiveness  of  the  boundary  condition.  Auxiliary  variables  (see  [Ref.  26])  will 
alleviate  this  problem. 

•  Distributing  the  C/s  on  the  interval  [.5, 1]  seemed  to  reduce  spurious  reflection, 
but  how  to  distribute  these  values  bestcouldnotbe  determined.  In  any  event, 
it  appeared  that  at  least  one  of  the  values  should  be  equal  to  the  gravity  wave 
speed  Co,  which  in  this  case  was  1. 

•  As  ^j— approached  5%,  spurious  reflections  became  visible.  Values  of  3% 

|^7 1  max 

produced  acceptable  results.  Obviously,  less  is  better. 

2.  Varying  U  and  V 

In  this  set  of  trials  U  and  V  are  varied,  however,  their  respective  values  are 
kept  the  same  in  each  layer.  The  following  parameters  are  fixed: 

Hi  =  {.03,  .07},  pi  =  {1.,  1.05}, 

At  =  .1,  Ax  =  Ay  =  .5,  (VII, 34) 

J=  3,  Cj  =  {.8,  .9, 1}. 

Results  are  posted  in  Table  VI.  These  trials  were  repeated  for  Ax  =  Ay  =  .25  with 
results  reported  in  Table  VII.  Another  set  of  trials  were  conducted  for  Cj  =  {.1,  .6, 1} 
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Table  VI.  Varying  U  and  V  with  Ax  =  Ay  =  .5  (U,  V  Equal  in  each  Layer) 


V 

V 

MU- 

1  T?l  rrma.3- 

{.0001, .0001} 

{-.0001, -.0001} 

1.92  x  10~a 

4.13% 

(.01, .01} 

{-.01, -.01} 

1.92  x  10-s 

4.13% 

(.025, .025} 

{-.025,-. 025} 

1.92  x  10-Q 

4.13% 

(.03,.  03} 

{-.03, -.03} 

1.92  x  10 -l3 

4.14% 

(.035,  .035} 

{-.035, -.035} 

1.92  x  10“° 

4.17% 

{.04, .04} 

{-.04, -.04} 

1.92  x  10"a 

4.20% 

{.05,.  05} 

{-.05, -.05} 

1.92  x  10-Q 

4.31% 

{•1,1} 

{-1,-1} 

1.93  x  10 -a 

7.01% 

{.25, .25} 

{-.25, -.25} 

2.32  x  10-3 

95.5% 

Table  VII.  Varying  V  and  V  with  Ax  =  Ay  =  .25  ( U ,  V  Equal  in  each  Layer) 


U 

V 

MU- 

ItZlrrmai 

{.0001, .0001} 

{-.0001, -.0001} 

1.42  x  10-° 

2.28% 

{.01, .01} 

{-.01,-01} 

1.42  x  10"a 

2.24% 

{■1-1} 

{-1,-1} 

1.40  x  10"6 

3.98% 

{.25, .25} 

{-.25, -.25} 

1.07  x  10-3 

58.5% 

with  Ax  =  Ay  =  .25.  All  other  parameters  were  unchanged.  Again  U  and  V  was 
varied,  but  this  time  the  advection  coefficients  were  allowed  to  differ  between  the 
layers.  Table  VIII  reports  the  results  when  the  disturbance  is  initiated  in  TL,  while 
Table  IX  reports  results  when  the  disturbance  is  initiated  in  L^.  From  these  results, 
it  can  be  concluded: 

•  The  problem  became  unstable  when  U,  V  became  large.  A  good  rule  of  thumb 
is  to  keep  U,  V  <£  d,  the  total  depth  of  the  medium.  This  is  a  physical  con¬ 
straint  in  an  ocean  environment  that  assumes  that  the  magnitude  of  advection 
will  be  much  less  than  the  magnitude  of  depth. 

•  Refining  the  grid  improved  performance. 

•  The  instability  of  the  problem  intensifies  when  U,  V  differs  between  layers, 
Such  problems  should  be  avoided,  when  advectdve  differences  between  layers 
are  small. 
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Table  VIII.  Varying  U  and  V  with  Ax  =  Ay  =  .25  for  Event  Initiated  in  Ll  (U,  V 
not  Necessarily  Equal  in  each  Layer) 


u 

V 

l^lfnoa! 

11 -II  mu 

Mn.« 

{.02, .02} 

{-.02, -.02} 

1.92  x  10-° 

3.50% 

{0,0} 

{-.02, -.02} 

1.92  x  10-“ 

3.51% 

{.02, .02} 

{0,0} 

1.92  x  10-° 

3.60% 

{0,0} 

{0,0} 

1.92  x  10“® 

3.92% 

{.01,0} 

{-01,0} 

1.92  x  10-® 

7.26% 

{0,-02} 

{0,-02} 

1.92  x  10“® 

12,6% 

{.01,- .01} 

{-.01,01} 

1.94  x  10-® 

12.8% 

{.02,0} 

{-.02,0} 

1.92  x  10“® 

13.1% 

{.02,- .02} 

{-.02, .02} 

1.99  x  10-® 

25.5% 

3-  Varying  Layer  Thicknesses  in  a  2-Layer  Problem 

In  this  set  of  trials,  the  layer  thicknesses  are  varied.  The  following  parameters 


were  fixed: 

Pi  =  {  1..1.05}, 

At  =  .1  Ax  =  Ay  =  .5 

J  =  3  Cj  =  {.8,  .9, 1} 

U  =  {.025,  .025},  V  =  {-.025,  -.025}. 


(VII. 35) 


Table  DC.  Varying  U  and  V  with  Ax  =  Ay  =  .25  for  Event  Initiated  in  L 2  (U,  V  not 
Necessarily  Equal  in  each  Layer) 


V 

V 

{0,0} 

{-.02, -.02} 

3.61% 

{01,-01} 

{-.01,01} 

3.63% 

{0,02} 

{0,-02} 

1.92  x  10-® 

3.71% 

{.02, .02} 

{-.02, -.02} 

1.92  x  10-® 

3.74% 

{0,0} 

{0,0} 

1.92  x  10-® 

4,38% 

{.02,- .02} 

{-.02, .02} 

■KiliCTliBa 

6.31% 

{.02,  .02} 

{0,0} 

6.32% 

{.01,0} 

{-.01,0} 

6,85% 

{02,0} 

{-02,0} 

1.92  x  10-® 

11.3% 
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Results  for  a  disturbance  initiated  in  Li  and  i2  are  posted  in  Table  X.  It  is  readily 

Table  X.  Varying  Layer  Thickness  in  a  2- Layer  Problem  with  an  Li  Initiated  Event 
(left)  and  L2  Initiated  Event  (right) 


concluded  that; 

•  Varying  layer  thicknesses  does  not  alter  the  effectiveness  of  the  Higdon  NRBC. 
This  conclusion  is  independent  of  the  layer  in  which  the  disturbance  was  ini¬ 
tiated, 

4.  Varying  Density  Distribution  in  a  2-Layer  Problem 

In  this  set  of  trials,  the  density  distributions  are  varied.  The  following  param¬ 
eters  were  fixed; 

Oi  =  {.03,  .07}, 

At  =  .1  Ax  =  A  y  =  .5 

(VII, 36) 

J  =  3  Cj  =  {.8,  .9, 1} 

U  =  {.025,  .025},  V  =  {-.025,  -  .025}. 

Results  for  a  disturbance  initiated  in  Li  and  L2  are  posted  in  Table  XI,  From  these 
trials  it  is  concluded  that; 

•  Density  changes  do,  to  some  extent  alter  the  behavior  of  the  Higdon  HNRBC, 
but  in  general,  they  remain  effective  as  long  as  the  density  increases  monoton- 
ically  with  each  increasing  layer. 
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Table  XI.  Varying  Layer  Thickness  in  a  2- Layer  Problem  with  a  Ll  Initiated  Event 
(left)  and  L2  Initiated  Event  (right) 


pi 

\fl\max 

||g||mair 

Mmn 

|?7  Imaa 

llgl|may 

{1,1} 

1.89  x  10-s 

3.58% 

1.89  x  10-s 

3.58% 

{1,1.05} 

1.92  x  10"a 

3.47% 

1.89  x  10-° 

3.58% 

O 

i—| 

i — 1 

i—| 

1.99  x  10-s 

3.20% 

1.88  x  10-s 

3.60% 

{1,1.15} 

2.03  x  10"a 

2.02% 

1.88  x  10"a 

3.60% 

{1,1.20} 

2.05  x  10~a 

2.98% 

1.88  x  10"s 

3.57% 

{1,1.25} 

2.05  x  10-° 

3.02% 

1.91  x  10"° 

3.47% 

{1,1.50} 

2.18  x  10-° 

3.11% 

2.10  x  10"° 

3.20% 

{1.05,1} 

1.85  x  10"s 

4.35% 

1.89  x  10~a 

3.92% 

{1.10,1} 

1.92  x  10-a 

8.84% 

1.93  x  10"a 

3.92% 

{1.25,1} 

2.10  x  10-6 

8.69% 

1.09  x  10-6 

6.94% 

5.  Varying  At  in  a  2- Layer  Problem 

In  this  set  of  trials  At  is  varied,  The  following  parameters  were  fixed: 


fit  =  {.03,  .07},  ^  =  {1,1.05} 
Ax  =  Ay  =  .5, 

J  =  3,  Cj  =  {.8,  .9, 1}, 

V={  0,0},  V  =  {0, 0}. 


(VII.37) 


Results  are  posted  in  Table  XII. 


In  general,  decreasing 


At  does  not  change  the 


Table  XII.  Varying  At  in  a  2- Laver  Problem  with  a  Initiated  Event 


At 

|^7 1  max 

MruM 

.1 

1.92  x  10^ 

4.13% 

.05 

3.82  x  10"° 

3.18% 

.01 

1.91  x  lO"6 

3.25% 

.005 

3.82  x  10-5 

3.23% 

.002 

9.57  x  10-5 

3.23% 

effectiveness  of  the  Higdon  NRBC.  However  it  appears  that  1^1™^  is  inversely  pro¬ 
portional  to  At.  The  probable  reason  for  this  behavior  is  that  the  introduction  of 
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the  event  produces  a  discontinuity  in  rj  with  regards  to  t.  A  smaller  At  results  in  a 
more  pronounced  discontinuity.  It  takes  several  numerical  time-steps  to  sort  out  the 
discontinuity.  The  result  is  a  lr/1^^  that  is  greater  then  the  maximum  amplitude  of 
the  event  induced  on  12. 

G.  THE  “HIGDON  MATRIX” 


r  ln^j«  tar  A  «ata  lUJ 


X  IDD  I3D  XD  an  XD  HD  -&□ 


Figure  60.  Higdon  Matrix  Image  for  £2  (20  x  20)  with  Higdon  NRBC’s  with  Order 
J  =  9  Applied  to  Four  Sides 

As  mentioned  in  Section  VILA,  when  non-zero  advection  terms  are  incorpo¬ 
rated  into  the  shallow  water  model  the  problem  must  be  solved  implicitly,  and  a 
N£  X  ATy  matrix  with  a  bandwidth  of  2A^  is  generated.  An  image  of  this  matrix 
is  presented  in  Figure  60  where  zero  elements  are  black  and  non-zero  elements  are 
white.  Here  the  truncated  domain  £2  is  approximated  using  21  x  21  grid  and  Higdon 
NRBC’s  of  order  J  =  9  are  applied  to  all  four  sides.  On  the  top  and  the  bottom  of 
the  image  we  see  10  light  diagonal  lines.  These  lines  represent  the  discretization  for 
the  y-boundaries  T /v  (top)  and  Tjv  (bottom).  The  heavier  line  along  the  diagonal  is 
three  points  thick  and  is  flanked  to  the  left  and  right  by  two  thinner  lines.  These 
result  from  the  discretization  of  the  interior  points.  Finally,  the  periodic  “short-spike” 
pointing  to  the  left  and  right  were  generated  by  the  Higdon  NRBC’s  on  and  rH,- 
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respectively.  Note  that  there  are  only  19  each  of  these  short  horizontal  lines,  This 
indicates  that  the  corner  points  were  included  in  the  ^-boundaries,  otherwise  21  (Nx) 
of  such  pairs  would  "be  visible. 

It  is  evident  from  the  image,  that  the  Higdon  matrix  required  for  non-zero 
advection  problem  is  sparse.  The  number  of  non-zero  points  generated  by  the  domain 
interior  is: 

5W,  -  2)(ATy  -  2),  (VII. 38) 

and  the  number  of  non-zero  points  generated  by  the  four  Higdon  NRBC’s  is: 

(2NX  -|-  2Ny  —  4)( J 1).  (VII. 39) 


Therefore  the  fraction  of  non-zero  elements  in  the  matrix  is: 


5(7Va  -  2)(7Vy  -  2)  +  (2NX  4-  2TVy  -  4 )( J  +  1) 


(VII  .40) 


In  the  case  of  our  example  where  Nx  =  Ny  =  21  and  J  =  9,  only  1.34%  of  the  matrix 
is  populated  with  non-zero  values.  In  the  TV-layer  model,  N  of  such  matrices  might 
be  produced  severely  taxing  computer  memory.  Increasing  domain  size  and  a  refining 
grid  would  further  exacerbate  the  problem.  Clearly  sparse  matrix  procedures  are  in 
order  and  should  be  investigated  to  stream  line  the  non-zero  advection  scheme  and 
produce  a  faster  algorithm. 
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VIII.  CONCLUSION  AND 
RECOMMENDATIONS  FOR  FURTHER 

RESEARCH 

From  the  preceeding  investigation  we  conclude  that  a  very  large  domain  with 
dispersive  wave  action  governed  by  linear  SWE’s  can  be  effectively  restricted  using 
Higdon  NRBC’s.  This  boundary  condition  will  also  work  for  versions  of  the  linear 
SWE  that  include  the  effects  of  stratification  and  advectdon.  However,  as  the  SWE 
model  becomes  more  complex,  careful  consideration  must  be  given  to  values  assigned 
to  problem  parameters  to  ensure  the  stability  of  the  scheme. 

There  are  many  aspects  of  the  problem  addressed  in  this  dissertation  that 
should  be  investigated  further.  With  regards  to  the  SWE  and  the  geophysical  envi¬ 
ronment,  further  model  development  should  include: 

•  Discretization  of  the  non-linear  SWE. 

•  Incorporation  of  a  non-constant  Coriolis  parameter. 

•  Inclusion  of  a  bottom  topography  that  varies  with  x  and  y. 

•  Consideration  of  terms  resulting  from  the  Earth’s  curvature  (i.e.  domains  with 
horizontal  dimensions  greater  than  1000  km). 

With  regards  to  computational  techniques,  one  should  consider  using  sparse  matrix 
algorithms  for  models  requiring  implicit  solution  techniques.  Finally,  with  respect  to 
the  Higdon  NRBC,  the  following  areas  of  research  are  recommended: 

•  Development  of  schemes  that  use  auxiliary  variables  to  eliminate  non-zero 
buffer  zones.  This  is  done  outside  this  dissertation. 

•  Development  of  schemes  to  optimize  selection  of  Higdon  coefficients. 

•  Development  of  methods  to  update  Higdon  coefficients  dynamic  ally/  adaptively. 

•  Extension  of  the  rectangular  domain  to  three  dimensions. 

•  Development  of  Higdon  NRBCs  for  cylindrical  and  spherical  coordinates. 
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•  Creation  of  two-way  Higdon  NRBC’s  to  incorporate  data  from  the  immediate 
vicinity  of  the  truncated  domain. 

These  extensions  to  the  current  research  are  potentially  of  great  value  to  oceanogra¬ 
phers  and  meteorologists  alike  and  may  further  promote  the  use  of  Higdon  NRBC’s 
in  weather  prediction  models. 
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